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CHAPTER  1 

INTRODUCTION 


Central  to  a  broad  class  of  problems  is  the  state  estimation  of  a 
hybrid  Markov  process,  for  which  the  underlying  process  consists  of  a 
Markov  chain  driving  continuous  dynamics.  Examples  wherein  this  issue 
arises  include  failure  detection  problems,  where  a  system  subject  to 
failures  can  be  considered  as  a  single  hybrid  Gauss-Markov  process,  and 
multiobject  tracking,  where  each  of  several  targets  car  be  modelled  as 
a  hybrid  Gauss-Markov  process. 

The  purpose  of  this  work  is  to  derive  algorithms  for  estimating 
the  state  sequence  of  single  and  multiple  hybrid  Markov  processes,  using 
all  available  measurement  data.  A  maximum  a  posteriori  (MAP)  decision 
criterion  is  used  for  the  choice  of  the  estimates,  which  leads  to  some 
nice  recursive  implementations.  In  the  case  of  linear,  Gaussian  dynamics, 
these  will  take  the  form  of  coupled  Viterbi  and  Kalman  algorithms. 

Since  limited  resources  are  available  for  the  implementation  of 
the  algorithms,  the  amount  of  computation  must  be  kept  at  a  reasonable 
level:  this  suggests  the  use  of  some  techniques  for  reducing  the  com¬ 
putational  burden.  While  many  ad  hoc  methods  for  doing  this  have  been 
proposed  (see  [9]-[10]),  this  thesis  derives  techniques  which  achieve 
this  goal  with  absolutely  no  increase  in  error  probability. 

For  this  purpose,  we  proceed  step  by  step.  Chapter  Two  will 
specify  the  mathematical  model  to  be  used,  and  Chapter  Three  the  re¬ 
lated  work  done  in  this  area.  In  Part  II,  the  single  system  problem 
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is  discussed:  Chapter  4  reviews  the  discrete  state  case  only.  Chapter  5 
the  continuous  state  case  only,  and  the  combination  for  a  hybrid  state 
estimation  is  developed  in  Chapters  6-7.  Then,  Fart  III  studies  the 
multiple  hybrid  system  case.  The  main  issue  in  this  part  is  to  obtain 
a  partial  decomposition  of  our  estimation  problem  into  several  indepen¬ 
dent  subproblems:  this  is  the  notion  of  clustering,  which  greatly  re¬ 
duces  the  computational  complexity.  An  application  to  multitarget 
tracking  is  discussed  in  Chapter  11,  where  we  prove  the  existence  of  an 
optimal  gating  leading  to  the  decomposition  of  a  multitarget  tracking 
problem  into  independent  single  target  tracking  problems. 
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CHAPTER  2 

FORMAL  PROBLEM  STATEMENT 


In  this  chapter,  we  formulate  the  MAP  estimation  problems  for  the 
single  and  multiple  hybrid  Markov  processes  in  relation  to  the  two  basic 
examples  that  motivate  them:  failure  detection  and  multitarget  tracking. 

Since  the  purpose  of  our  work  is  to  derive  computational  algorithms 
all  events  are  assumed  to  occur  at  discrete  instants  of  time  k  »  0,1,2... 
Throughout  this  work,  the  use  of  the  superscript  (  )  will  denote  the 

k 

sequence  up  to  and  including  k:  a  ■  (a(0) , . . . ,a(k)) . 

2.1  Necessity  for  a  Hybrid  State  Description 

In  many  estimation  problems,  the  description  of  the  state  of  a 
dynamic  system  only  in  terms  of  a  continuous  state  is  not  sufficient. 

For  instance,  in  failure  detection  problems,  the  system  considered 
is  subject  to  abrupt  changes  (failures  or  repairs).  The  normal  opera¬ 
tion,  or  no  failure,  model  of  the  system  is  described  in  standard  state 
space  form: 

system  dynamics:  x(k+l)  *  F  x(k)  +  G  u(k)  +  wQc) 
sensor  equation:  t(k)  ■  H  x(k)  ♦  v(k) 

The  abrupt  changes  can  arise  in  a  number  of  ways;  they  can  manifest 
themselves  as  actuator  failures:  shift  in  the  control  gain  matrix  G  or 
increased  process  noise;  sensor  failures:  abrupt  changes  in  H  or  increase 
in  measurement  noise;  or  as  actual  changes  in  the  plant  dynamics  F.  The 
mode  (failed  or  unfailed)  under  which  the  system  operates  at  a  given 
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time  is  a  discrete  state  complementing  the  usual  continuous  state,  whose 
dynamics  can  be  described  by  a  state-transition  diagram  (see  Fig.  2.1). 
This  combined  hybrid  state  provides  the  most  complete  description  of 
the  system. 

Similarly,  in  multitarget  tracking  problems,  it  seems  natural  to 
describe  the  set  of  targets  with  a  hybrid  process.  In  an  area  under 
surveillance,  an  unknown  number  of  targets  are  moving;  some  new  targets 
can  arrive  suddenly,  while  others  can  disappear.  This  process  can  be 
simply  characterized  using  an  augmented  state  for  each  target.  We  as¬ 
sume  a  large  but  fixed  number  of  targets  in  the  area,  many  of  which  are 
dead  or  in  discrete  state  0.  A  new  arrival  is  the  transition  for  a 
target  between  state  0  and  the  alive  discrete  state  1;  a  disappearance 
is  the  reverse  transition.  Thus,  each  target  is  described  with  a  hybrid 
state:  discrete  state  0  or  1  and  usual  continuous  states;  e.g.,  positions 
and  velocities.  With  more  complex  discrete  state  models,  a  variety  of 
phenomena  such  as  maneuvers ,  signal  fading,  and  variable  configurations 
can  be  modeled  (see  Fig.  2.2). 

In  both  examples,  we  see  that  the  introduction  of  a  hybrid  state 
description  for  the  system  provides  the  most  compact  and  complete  charac¬ 
terization  of  the  process. 

We  next  formulate  the  single  hybrid  state  estimation  problem. 

2.2  The  Single  Hybrid  State  Estimation  Problem 

In  this  section,  we  develop  the  single  hybrid  Markov  process  model 
and  formulate  our  objective. 


T(juo  :  4  Mocnefe  sfoik> 


00 :  ackuj-ac.  <aW  jb&KDOt  ma^cu-SscL 
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The  state  s(k)  of  a  single  hybrid  Markov  process  at  time  k  belongs 
to  a  state  space  S.  It  consists  of  both  a  discrete  state  d(k)  and  a 
continuous  state  x(k) .  The  discrete  state  d(k)  belongs  to  a  finite  set 
D  ■  (l, . . . ,n) .  The  continuous  state  xfk)  belongs  to  a  state  space  X  of 
finite  dimension  n. 

The  dynamics  of  s(k)  are  assumed  to  follow  a  Markov  process,  in 
the  sense  that  the  probability  distribution  P(s (k+1) !s(k) , . . . ,s(0))  of 
being  in  state  s(k+l)  at  time  k+1,  given  all  states  up  to  time  k,  de¬ 
pends  only  on  the  state  s(k)  at  time  k: 

P(s(k*l)|sk)  -  P(sOk*l)|sCk))  (2.1) 

Furthermore,  the  discrete  state  dCk+1)  is  assumed  to  depend  only  on  the 
discrete  state  d(k)  at  time  k.  So,  P(s (Tc+1)  |  s(k]}  can  be  written  as: 

P(s(k+l)|s(k))  =  P(x(k+l)|x(k),d(k),d(k*l))*  P(d(k*l)|d(k)) 

‘(2.2) 

The  system  is  observed  during  an  interval  [0,K],  and,  for  simplicity, 
we  assume  that  the  measurement  process  is  continuous,  although  all  the 
results  can  be  generalized  to  a  hybrid  measurement  process  at  the  cost 
of  more  cumbersome  notation. 

At  each  discrete  instant  of  time  k,  a  continuous  measurement  z(k) 
of  the  hybrid  state  s(k)  is  generated.  The  observation  z(k)  belongs 
to  a  state  space  Z  of  finite  dimension  m.  The  measurement  process  is 
assumed  to  be  memory less:  the  probability  distribution  for  the  measure¬ 
ment  z(k)  depends  only  on  the  state  of  the  system  at  time  k: 

P(z(k)|sk,zk'1)  -  P(2(k) Is (k))  (2.3) 
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A  further  assumption  often  made  on  dynamic  systems,  when  the 
system  equations  can  be  linearised  and  a  good  approximation  of  the 
physical  noises  are  given  by  Gaussian  densities  in  the  Linear- Gaussian 
assumption.  Under  this  assumption,  the  distributions  P(x(k*l) | x(k) , 
d(k)  *  p,d(k+l)  *  q)  and  P(z(k)  |  x(k)  ,d(k)  *  p)  are  defined  by  the  fol¬ 
lowing  relations: 

x(k+l)  «  F(p,q,k)xCk)  *  w(p,q,k)  (2.4) 

z(k)  »  H(p,k)x(k)  +  v(p,k)  (2.S) 

F(p,q,k),  H(p,k)  are  known  n  x  n  and  ax  n  time-varying  matrices,  that 
depend  also  on  the  discrete  state  of  the  system.  w(p,q,k),  v(p,k) , 
x(d(0)]  are  independent  zero-mean  White  Gaussian  processes  with: 

E{w(p,q,k)wT(p,q,k)}  *  Q(p,q,k)  >  0 

E(v(p,q,k)vT(p,q,k)}  -  R(p,q,k)  >  0  (2.6) 

Efx[d(0)]xT[d(0)]}  =  P(d(0) )  >  0 

In  a  failure  detection  application,  the  equations  under  the  normal 
mode  (d(k)  *  0)  axe: 

x(k+l)  -  F  x(k)  +  w(.k) 
z(k)  -  H  x(k)  +  Vq (k) 

Then,  if  a  sensor  failure  arises  at  time  k  (d(k)  *  1),  the  measurement 
equation  becomes: 

zOO  ■  v^k) 

Therefore:  H(0,k)  «  H,  vC0,k)  -  vQ(k),  H(l,k)  ■  0,  v(l,k)  -  v^k). 
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The  model  developed  above  is  the  basis  for  Part  II  of  our  work, 
which  discusses  the  problem  of  estimating  the  state  of  a  single  hybrid 
Markov  process.  There  can  be  two  possible  objectives  in  this  estimation, 
depending  on  the  application. 

For  instance,  in  failure  detection  problems,  if  the  system  has  suf¬ 
ficient  hardware  back-up  capabilities,  then  the  objective  of  the  failure 
detection  system  is  to  quickly  detect  the  failures:  this  is  done  by  esti¬ 
mating  the  discrete  state  sequence  of  the  system.  But,  if  there  is  not 
enough  redundancy  in  the  system,  then  the  system  has  to  operate  under 
the  degraded  mode.  In  this  case,  the  objective  of  our  failure  detection 
system  is,  at  the  same  time,  to  detect  the  failures  and  to  perform  an 
estimation  of  the  continuous  state. 

So,  one  possible  objective  is  the  MAP  estimation  of  the  discrete 
It  K 

state  sequence  d  ,  based  on  z  ,  that  minimizes  the  probability  of  an 

K 

error  in  estimating  d  .  The  problem  is  to  find  d  such  that: 

P(dK|zK)  >  P(dK|zKj,  for  all  dK  G  DK  (2.7) 

In  other  cases,  where  we  are  interested  in  the  estimation  of  the  hybrid 

K  K 

state  sequence  s  ,  the  objective  is  the  MAP  estimation  of  s  ,  based  on 

zK.  The  problem  is  to  find  S*'  such  that: 

PC^jz*)  >  P(sK|zK),  for  all  sK  G  SK  (2.8) 

In  both  cases,  an  algorithm  for  estimating  the  state  sequence  (dis¬ 
crete  or  hybrid)  of  the  single  hybrid  system  must  keep  the  amount  of 
computation  at  a  reasonable  level,  since  limited  computational  resources 
are  available. 


We  next  consider  the  more  complex  multiple  hybrid  state  estimation 
problem. 

2.3  The  Multiple  Hybrid  State  Estimation  Problem 

In  this  section,  we  present  the  multiple  hybrid  Markov  process 
model,  formulate  our  objective,  and  see  how  it  relates  to  the  multi¬ 
target  tracking  problem. 

The  multiple  hybrid  system  consists  of  a  finite  collection  of 
single  hybrid  Markov  processes.  £■  *•  state  s(k)  at  time  k  is  equal  to: 

s(k)  •  (s, (k),..  w^qre  for  all  i  *  l,...,p, 

i  5* 

the  state  s^(k)  *  (d^(k) ,/ .  rk'’'f  :-i  subsystem  i  belongs  to  S. 

The  subsystems  are  is:  vied  to  have  independent  dynamics.  There¬ 
fore,  the  probability  distribution  P(s£k*l)|  s(k))  is  equal  to: 

.  P 

P(s(k+1)  |s(k))  *  H  P(s.Ck+l)!s.Ck))  (2.9) 

i»l  1 

and 

PCs.fk+Dls.Ckj)  »  PCxiCk+l)idi(k),d.Ck*l),x.Ck)j  x 

P(di(k+l)|diCk))  (2.10) 

As  mentioned  in  Section  1,  in  multitarget  tracking,  the  discrete 
state  of  a  target  d^k)  may  take  the  value  1  if  the  target  is  present 
in  the  area  at  time  k,  and  0  if  it  is  absent.  The  continuous  state 
x^k)  consists  of  the  positions  and  velocities  of  target  i.  The  inde¬ 
pendent  dynamics  assumption  means  that  we  consider  independent  positions 
and  velocities  of  targets  and  independent  arrivals  or  disappearances. 
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The  independent  arrival  assumption  is  certainly  the  most  difficult  to 
accept,  since  an  arrival  may  be  less  likely  to  occur,  if  there  are  al¬ 
ready  500  targets  in  the  area,  rather  than  only  5  targets. 

The  continuous  measurement  generation  process  for  the  multiple 
hybrid  system  is  the  following:  each  subsystem  i  generates  a  set  of 
measurements  Y.  (k)  at  time  k:  if  Y.  (k)  ■  {  z.  ?(k),...,z.  Ck)},  then 

1  1  X  ,  i  X 

for  all  j  ■  l,...,q,  z.  .  (k)  belongs  to  Z.  We  can  have  multiple  detec- 
tions  corresponding  to  the  same  subsystem  in  the  case  of  multiple, 
separate  sensors.  We  allow  Y^(k)  to  be  the  empty  set,  in  which  case  we 
say  that  subsystem  i  has  not  been  detected  at  time  k.  Besides,  because 
of  the  noisy  environment,  some  measurement  data  may  appear  independent 
from  all  of  the  system  states:  this  set  of  measurements  is  the  set  of 
false  alarms  and  is  denoted  by  Y^(k) . 

In  a  multitarget  environment,  a  missed  detection  can  arise  because 
of  a  sensor  failure,  or  some  natural  obstacle  between  a  sensor  and  the 
target.  A  false  alarm  could  come,  for  instance,  from  a  flock  of  geese, 
when  the  surveillance  system  is  observing  aircrafts. 

Each  subsystem  is  assumed  to  generate  its  set  of  measurements  inde¬ 
pendently  of  the  other  subsystems  and  of  the  false  alarms.  This  is  a 
reasonable  assumption  in  multitarget  tracking.  Therefore: 

P  (Y  ,(k) , . . . ,  Y  (k) ,  Y-Ck)ls(k))  -  n  P(Y.(k)[s.Ck))*P(YfCk)) 

i  pi  ii  i 

(2.11) 

where  for  all  i  *  l,2,...,p,  P(Y^(k) | s^(k))  and  PCYjCk))  are  assumed 
to  be  known. 

As  for  the  single  system  case,  we  have  an  additional  Linear-Gaussian 
assumption.  For  simplicity,  we  assume  that  only  one  sensor  is  used.  In 


this  case,  each  subsystem  i  generates  at  most  one  measurement  z^(k)  at 
time  k.  For  each  subsystem  i  =  l,...,p,  PCx^k+l) \  (k) ,  d^(k)  =  q, 

d^(k+l)  *  r)  and  P(z^(k) | x^(k) ,  d^(k)  3  q)  are  defined  as  for  the  single 
system  model  by: 


x^+1)  ■  F-^q.r.k^Ck)  +  w^k)  (2.12) 

zi(k)  «  H.fq.kJXj^Ck)  ♦  v.,(k)  (2.13) 

with 

E( w^(q,r,k)wT(q,r,k)}  *  Q^q.r.k) 

E(vi(q,r,k)vT(q,r,k)}  ■  R-^q.r.k)  (2.14) 

E(x.[d(0)]xj[d(0)]}  -  P.(d(0)) 


The  false  alarm  process  is  such  that  each  false  alarm  z^(k)  is  genercted 
independently  of  the  other  false  alarms;  we  will  see  in  part  III,  why 
we  do  not  assume  a  Poisson  distribution  for  the  number  of  false  alarms. 
The  probability  distribution  for  a  false  alarm  z^(k)  is  a  zero-mean 
Gaussian  density  with  a  known  covariance  Z . 

For  the  multitarget  tracking  application,  if  we  assume  planar 
straight-line  motions  and  position  measurements,  then  we  have: 


Xj(k) 


x\w~ 

vj(k) 
x?(k) 
v?(k)  _ 


1  2 

in  the  plane  (x  ,x  ), 


and 
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xi(k+l)  = 


^00  - 


1 

0 

0 

0 

1 
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1 

1 

0 

0 

0 

0 


0 

0 

1 

0 

0 

1 


0 

0 

1 

1 

0 

0 


x.Ck)  *  w.(k) 


x.Ck)  +  v.(k) 


The  problem  of  estimating  the  discrete  or  hybrid  state  sequence  of 
the  multiple  system  is  much  more  complex  than  the  simple  juxtaposition 
of  p  estimation  problems,  one  for  each  subsystem. 

In  the  multitarget  case,  the  objective  of  the  surveillance  system 
is  track  the  existing  targets,  to  announce  the  new  targets  and  the  dis¬ 
appearances  of  old  targets.  But,  because  of  the  cluttered  environment 
the  tracker  does  not  know  which  measurement  originated  from  a  given  tar¬ 
get  and  which  ones  are  the  false  alarms.  Therefore,  there  exists  an 
uncertainty  in  the  origin  of  the  measurements ;  and  any  target  state 
estimation  has  to  take  into  account  this  uncertainty. 

For  the  general  estimation  problem,  the  actual  set  of  observations 
is  denoted  by  Y(k)  *  { (k) , . . . ,  z^Ck)} .  The  uncertainty  is  introduced 
by  a  physical  permutation  P(k),  which  is  a  map  from  the  multiple  system 
and  the  set  of  false  alarms  onto  the  elements  of  Y(k).  This  permutation 
is  unknown;  therefore,  a  hypothesis  on  it  is  made  by  the  estimator,  that 
is  called  a  data  hypothesis  at  time  k  and  is  denoted  by  ACk) .  Each  data 
hypothesis  A(k)  has  the  following  features: 

-  A(k)  is  a  mapping  from  the  multiple  system  sCk)  into  the  set 


of  observations  Y(k), 


2' 


w 


-  Each  subsystem  i  is  associated  with  zero,  one,  or  more  measure¬ 
ments  of  Y(k).  A  subsystem  associated  with  no  measurement  is  naturally 
called  a  missed  detection  under  A(k) . 

-  Each  measurement  of  Y(k)  can  come  from  zero,  one,  or  more  sub¬ 
systems  of  s(k).  A  measurement  associated  with  no  subsystem  is  consi¬ 
dered  as  a  false  alarm  under  A(k) .  We  illustrate  this  mapping  in  Fig. 
2.3. 

We  assume  that  all  the  associations  of  a  measurement  to  a  subsystem 
are  a  priori  likely.  Therefore,  the  probability  P[A(k)|  s(k))  of  a  data 
hypothesis  A(k)  at  time  k,  depends  only  on  the  numbers  of  detections, 
missed  detections,  and  false  alarms  under  A(k) .  Furthermore,  we  assume 
for  simplicity  that  P(ACk)!s(k)]  depends  only  on  the  discrete  state 
d(k),  i.e.  P(A(k)[s(k))  ■  PCA(k) [ d(k)) .  In  multitarget  tracking,  this 
means  that  the  probability  of  detecting  a  target,  and  the  probability 
of  a  false  alarm  are  uniform  in  the  area.  If  we  make  this  assumption, 
then  P(ACk)|d(k))  takes  the  form: 

PCAOO|d(k))  *  PrfM^  detections  under  ACk)}  *  PrfM^  missed 
detections  under  A(k)}  x  PrCM^  false  alarms  under  A(k)}. 

This  data  association  uncertainty  complicates  the  estimation  prob¬ 
lems  as  follows. 

As  for  the  single  system  model,  there  are  two  different  objectives 
in  a  multiple  hybrid  state  estimation  problem. 

The  first  possible  objective  is  the  MAP  estimation  of  the  hybrid 

!C  K  K  15  K 

state  sequence  s  ,  based  on  YN.  The  problem  is  to  find  s  £  (S*j  for 

K  JC  K 

which  P(s  |Yj,  or  equivalently,  II  P(YCk)|s(k))  x  PCs(k)  ]  s(k-l))  is 

k«0 
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maximum.  But  now,  P(Y(k)js(k))  must  be  evaluated  from: 


P(Y(k) | s(k) )  *  £  P(Y(k)|s(k),  A(k))  x  P(A(k) | d(k) ) ,  (2. IS) 

(A(k)} 

and 

,  P 

P(Y(k) | s (k) ,A(k))  *  JI  P(Y  (A(k))|s.(k))  x  P(Yf(A(k))),  (2.16) 

i«l 

where  the  summation  is  over  all  possible  data  hypotheses,  and 
(Y^AQO) , . . .  ,Yp(A(k))  ,Yj(A(k))}  is  the  hypothesized  partition  of  Y(k) 
under  A(k) . 

There  exists  also  a  second  possible  objective.  For  instance,  in 

some  multitarget  tracking  problems,  we  can  be  interested  in  first 

deciding  upon  the  presence  of  targets  in  the  area,  then  in  locating 

them  based  on  the  above  decisions.  In  this  case,  the  objective  is 

K  K 

rather  the  MAP  estimation  of  d  based  on  Y  ,  which  minimizes  the  proba- 

K  p  y 

bility  of  an  error  in  estimating  d  .  The  problem  is  to  find  d  G  (D- ) 

K 

for  which  P(dK|YK),  or  equivalently,  H  PCY(k)[dk,Yk  1)  x  P(d(k)  I  if*.  *1)) 

nz-vri.'*  1  jk  vk-l)  must  be  IvSluated  from: 
zs  maximum.  P(Y(k)jd  ,Y  1 

P(Y(k)|dk,Yk_1)  -  £  P (Y (k)  [  dk , Yk” 1 , Ak)  x  P(Akjdk)  (2.17) 

{  AK} 

P 

n 

i-l 

P(Yf(A(k))  .  (2.18) 


and 


P (Y (k) | dk ,Yk~ 1 , Ak)  -  n  P(Y. (A(k) ) [ d^ ,Yk~* (Ak”X) )  x 


Under  both  objectives,  we  must  find  an  algorithm  for  estimating 
the  state  of  the  multiple  system,  that  keeps  the  amount  of  computation 
at  a  reasonable  level,  with  respect  to  our  computational  resources. 
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2.4  Summary 

In  this  chapter,  we  have  developed  the  basic  models  and  objectives 
for  the  single  and  multiple  hybrid  Markov  processes,  in  relation  to  the 
failure  detection  and  multitarget  tracking  problems. 

We  next  review  some  previous  work  done  in  these  areas,  and  compare 
it  with  the  approach  taken  here. 
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CHAPTER  3 

BACKGROUND  AND  LITERATURE  SURVEY 


In  this  Chapter,  we  compare  the  formulation  and  objectives  of  the 
last  chapter,  with  some  previous  work  done  in  failure  detection  and 
multitarget  tracking. 

3. 1  Failure  Detection 

A  survey  of  design  methods  for  failure  detection  systems  has  been 
done  by  Willsky  in  [5]. 

One  of  these  methods  is  closely  related  to  our  approach:  this  is 
the  multiple  hypothesis  filter-detectors,  in  which  a  bank  of  Kalman 
filters  is  constructed,  based  on  the  different  hypotheses  concerning 
the  system  behavior.  This  method  computes  the  probability  of  each  type 
of  failure  under  consideration.  The  problem  is  the  exponentially 
growing  bank  of  filters,  so  that  some  ad-hoc  pruning  techniques  have  to 
be  used,  to  limit  the  computational  complexity:  unlikely  hypothesis  se¬ 
quences  are  simply  dropped  during  their  formation,  or  a  single  shift  in 
the  system  behavior  is  assumed  every  K  steps. 

The  problem  of  detecting  a  single  switch  between  two  dynamic  models 
(failed,  unfailed)  is  also  treated  in  [6J-[7],  using  a  sequential  proba¬ 
bility  ratio  test  (SPRT) .  The  SPRT  optimizes  a  combination  of  the  time 
to  reach  a  decision  and  the  probabilities  of  making  wrong  decisions. 

Thus  many  "optimal"  algorithms  for  detecting  failures  use  some  ad- 
hoc  procedures  in  order  to  meet  some  reasonable  time  and  storage  require- 
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ments,  at  the  expense  of  a  degradation  in  the  performance. 

We  have  seen  in  Chapter  2  that  the  single  hybrid  state  estimation 
problem  applies  to  the  failure  detection  problem.  In  particular,  a 
hypothesis  in  the  multiple  hypothesis  filter-detectors  corresponds  to 
a  discrete  state  sequence  in  our  formulation.  The  objective  of  part  II 
is  to  find  an  algorithm,  that  uses  truly  optimal  pruning  techniques  to 
reduce  the  necessary  amount  of  computation.  If  this  algorithm  meets  also 
the  storage  requirements,  then  we  would  get  an  optimal  algorithm  for  the 
failure  detection  problem. 

3.2  Multitarget  Tracking 

An  excellent  survey  of  recent  work  in  multitarget  tracking  methods 
(up  to  1978)  has  been  done  by  Bar-Shalom  [8] . 

An  important  paper  in  this  field  has  been  given  by  Reid  [9].  A  set 
of  data-association  hypotheses  that  take  into  account  all  possible  ori¬ 
gins  of  every  measurements  (prior  targets,  new  targets,  and  false  alarms) 
is  generated  and  organized  sequentially  into  a  tree.  For  each  hypothesis, 
a  set  of  Kalman  Filters  is  constructed  to  track  the  targets  as  implied 
by  that  hypothesis.  The  probability  of  each  hypothesis  is  computed  se¬ 
quentially,  using  the  predicted  value  of  the  measurement  obtained  from 
the  state  estimates  of  the  filters:  the  closer  a  measurement  is  to  a 
particular  prediction,  the  higher  the  probability  that  the  measurement 
is  associted  with  the  corresponding  hypothesis.  A  key  idea  to  reduce 
computation  is  clustering .  In  the  hypothesis  tree,  any  branch  whose 
probability  lies  below  a  threshold  is  dropped.  The  probability  of  a 
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hypothesis  below  the  threshold  is  equivalent  to  the  present  measurement 
falling  outside  a  gate  drawn  around  the  prediction.  If  the  gates  of  the 
various  data  hypotheses  do  not  overlap,  then  the  set  of  received  measure¬ 
ments  can  be  partitioned  into  several  subsets,  each  subset  containing  the 
measurements  falling  outside  a  gate;  and  the  hypothesis  tree  can  be  broken 
up  into  smaller,  independent  trees.  If  some  of  the  gates  do  overlap,  then 
a  tree  can  be  formed  for  the  cluster  of  the  corresponding  data  hypotheses. 
This  clustering  leads  to  a  great  simplification,  since  the  total  number 
of  hypotheses  in  the  smaller  trees  will  be  much  less  than  the  number  of 
hypotheses  in  the  global  tree.  A  new  target  is  hypothesized,  as  each  new 
measurement  is  received;  it  is  confirmed,  if  after  dropping  the  unlikely 
hypotheses,  the  target  has  a  probability  of  existing  equal  to  one.  The 
possibility  of  a  target  ceasing  to  exist  is  not  covered  by  the  paper. 

Thus,  this  algorithm  uses  some  suboptiaal  pruning  techniques,  in 
order  to  reduce  the  amount  of  computation,  and  some  ad-hoc  procedures 
to  decide  upon  the  arrival  of  a  new  target. 

Keverian  modifies  slightly  this  algorithm  in  [10]  and  Hughes  [11] 
applies  it  to  a  distributed  system. 

An  early  work  by  Sittler  [12J  proposed  to  compute  sequentially  a 
score  function,  using  the  inverse  of  a  likelihood  function;  a  new  track 
is  initiated  when  the  score  function  stays  above  a  threshold,  and  is 
doppred  when  it  goes  below  another  one.  A  modernized  version  of  this 
algorithm  has  been  done  by  Stein  and  Blackman  in  [13]. 

Morefield  {14]  proposed  a  solution  for  the  most  likely  data  associ¬ 
ation  in  an  integer-progTamming  framework;  it  is  a  batch  processing 
technique. 


27 


Alspach  [15]  computes  the  a  posteriori  density  of  the  state  of  the 
target,  which  is  a  Gaussian  Sum  density  because  of  the  different  pos¬ 
sible  data  hypotheses  to  this  target.  After  dropping  the  terms  in  the 
sum  with  small  weighting  coefficients,  the  state  estimate  is  found  as 
a  convex  combination  of  the  mean  of  the  individual  terms. 

Another  approach  is  taken  by  Bax-Shalom  [16]- [17]:  given  a  target 
state  has  been  initialized,  the  probability  that  a  current  measurement 
comes  from  that  target  is  computed.  Then  all  the  measurements ,  weighted 
by  these  probabilities,  are  used  to  update  the  target  state  estimate. 
This  is  the  Probabilistic  Data  Association  Filter.  No  initialization 
procedure  is  included  in  this  algorithm. 

We  have  seen  in  Chapter  2  that  the  multiple  hybrid  state  estima¬ 
tion  problem  applies  to  the  multitarget  tracking  problem.  In  our  for¬ 
mulation,  we  take  a  target-oriented  approach:  each  target  has  an  a 
priori  probability  of  appearing  and  disappearing,  then  the  received 
measurements  will  confirm  or  deny  its  presence  in  the  area.  Contrary 
to  Reid's  algorithm,  this  allows  a  systematic  procedure  for  track 
initiations  and  terminations,  since  the  MAP  hybrid  state  estimation 
algorithm  decides  upon  the  discrete  state  of  each  target.  To  keep 
the  amount  of  computation  at  a  reasonable  level,  this  algorithm  must 
use  some  pruning  techniques  in  the  tree  of  the  possible  discrete  state 
sequences  of  the  multiple  system.  We  are  interested  only  in  pruning 
rules  which  never  degrade  performance,  not  in  some  ad-hoc  heuristics; 
this  applies  to  part  II  as  well  as  part  III. 

Thus,  the  purpose  of  our  work  in  part  III  is  to  find  some  optimal 
pruning  rules  during  the  formation  of  this  tree,  and  particularly  to 
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examine  the  existence  of  optimal  gating,  that  could  greatly  reduce  the 
computational  complexity,  without  increasing  the  probability  of  an  error 
in  the  estimation. 

The  resulting  algorithm  for  estimating  the  multiple  state  sequence, 
even  if  it  exceeds  the  computational  requirements  using  these  rules  and 
gating  would  still  provide  a  basis  for  rational  (if  not  optimal)  proce¬ 
dures  for  deciding  upon  new  arrivals  and  disappearances  in  the  area,  and 
for  locating  the  targets. 

3.3  Summary 

In  this  chapter,  we  have  reviewed  some  previous  work  done  in 
failure  detection  and  multitarget  tracking,  and  we  have  compared  it 
with  our  formulation  and  objectives. 

We  next  consider  our  first  problem:  the  single  hybrid  state  esti¬ 


mation. 
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PART  II 


THE  SINGLE  SYSTEM  CASE 
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CHAPTER  4 

STATE  ESTIWTION  OF  A  MARXOV  CHAIN: 

THE  VITERBI  ALGORITHM 

In  this  chaptex,  we  describe  the  Viterbi  Algorithm  [1],  which  is 
a  recursive,  optimal  solution  to  the  problem  of  estimating  the  state  se¬ 
quence  of  a  finite  Markov  chain  observed  in  memoryless  noise.  This  algo¬ 
rithm  has  many  applications  to  digital  estimation  problems:  convolutional 
coding,  where  an  encoded  sequence  is  sent  through  a  memoryless  channel; 
digital  transmission  through  analog  channel. 

4.1  Problem  Statement 

In  this  part,  a  single  hybrid  state  system  is  considered.  The 
system  model  is  the  one  developed  in  Chapter  2,  Section  2.  Here,  the 
continuous  state  space  X  reduces  to  0,  so  that  only  the  discrete  state 
part  is  considered. 

The  Markov  and  memoryless  properties  (2.1)-(2.3)  become: 


P(d(k*l)ldk)  -  P(d(k+l)]d(k)) 

(4.1) 

P(z(k)|dk,zk“l)  -  P(z(k)|d(k)) 

(4.2) 

These  probabilities  may  be  time-varying,  but  we  do  not  explicitly  indi¬ 
cate  this  in  the  notation. 

The  system  is  observed  during  an  interval  [0,K], 

The  objective  is  to  minimize  the  probability  of  an  error  in  esti 

V 

mating  the  whole  sequence  d  »  (d(0) , . . . ,d(K)) •  The  solution  to  this 

problem  is  given  by  the  state  sequence  d  for  which  the  a  posteriori 
Ki  1C 

probability  P(d  j  is  maximum.  In  short,  the  problem  is  to  find  d 
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such  that: 

VdK  €  DK,  P(dK|zK)  >  P(dXjziS  . 

aV 

The  straightforward  computation  of  d  is  c:  *  Lnatorial :  find  the 

X 

most  likely  state  sequence  among  N  state  sequences. 

Using  the  Markov  and  memoryless  properties  (4.1)-(4.2),  we  can 
show  that  this  estimation  problem  is  formally  identical  to  the  problem 
of  finding  the  shortest  path  through  a  certain  graph,  which  requires  much 
less  computation.  This  is  the  basic  idea  in  the  Viterbi  Algorithm. 


The  Algorithm 


A  graph  can  be  associated  with  our  Markov  chain.  Each  node  cor- 
responds  to  a  distinct  state  at  a  given  time  and  each  branch  represents 
a  transition  to  some  new  state  at  the  next  instant  of  time.  The  graph 
begins  at  time  0  'with  a  state  dQ  that  is  assumed  to  be  known.  A  trivial 
extension  can  be  made  to  the  case  where  d(0)  is  given  by  some  a  priori 
probability  distribution  function.  To  every  possible  discrete  state 
sequence  d  ,  there  corresponds  a  unique  path  through  the  graph,  and 
vice-versa. 

K  K  K* 

The  problem  of  finding  d  for  which  P(d| z  j  is  maximum,  or  equi¬ 
valently  for  which  P(dX, zX)  =  PCd^jz^)  x  PCz^)  is  maximum,  is  the  same 
as  the  one  of  finding  the  path  in  the  graph  whose  length  XCd1^  * 

-in  P(dX,z^  is  minimum. 

Due  to  the  Markov  and  memoryless  properties  (4.1)-(4.2),  we  have: 

X 

P(dK,zK)  *  n  P(d(k)|dCk-l))  x  PCz(k)jdCk))  . 
k*0 

If  we  assign  each  branch  the  length  A(d(k) ,d(k-l))  ^  -in  P(d(k) j dfk-1) ) 
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-  in  P(z(k) | d(k)) ,  that  depends  on  the  received  measurement  z(k)  at  time 
k,  then: 


ACd*) 


K 

E  ACdCk).d(k-l)) 
k»0 


(4.3) 


the  total  length  of  a  path  is  equal  to  the  sum  of  the  length  of  the 


branches  composing  the  path. 

k 

In  the  graph,  d  corresponds  to  a  path  starting  at  the  node  dQ 
and  terminating  at  d(k) .  For  any  particular  node  d(k) ,  there  will  be 


several  such  paths.  The  shortest  path  from  d^  to  d(k)  is- called  the 
survivor  corresponding  to  d(k)  and  is  denoted  d[d(k)].  For  any  time 


k  >  0,  there  are  N  survivors  in  all,  one  for  each  dfk). 

-s 

The  shortest  complete  path  d  must  begin  with  one  of  these  sur- 

lc  k 

vivors.  This  is  true  because  of  the  following  property:  let  d^  and  d, 
be  two  state  sequences  up  to  time  k,  such  that  dj(k)  *  d2(k).  If 
A(d^)  <  X(d2),  then,  for  all  z(k+l)  , . . .  ,200  *  for  all  d(k+l) , . . .  ,d(K) , 
A(dJ)  <_  X (d^)  ,  where  d*  -  (d*,d(k*l) , . . . ,d(K))  and  d£  -  (d^dCk+l) , . . . , 
d(K)),  because  of  property  (4.3).  We  illustrate  this  property  in  Figure 


(4.1). 

R  a  {£ 

Therefore,  d2  cannot  be  a  part  of  the  shortest  path  d  and  can 
be  immediately  discarded  at  time  k. 


Thus,  only  the  N  survivors  dfd(k)]  and  their  length  T[d(k)]  * 
X[d(d(k))]  need  to  be  stored  at  time  k. 

The  recursion  proceeds  as  follows:  each  survivor  d[d(k)]  is  ex¬ 
tended  by  one  time  unit  and  the  length  of  the  extended  paths  are  com¬ 
puted,  using  the  received  measurement  at  time  k+1.  Then,  for  each  node 
d(k+l),  the  shortest  path  terminating  at  d(k+l)  is  selected  as  the  cor- 


o->  as  o-> 
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responding  survivor  d[d(k+l)].  This  procedure  is  repeated  until  k  *  K 
and  the  number  of  survivors  never  exceeds  N.  Thus,  at  time  K,  we  have 

A 

N  survivors  d[d(K)].  The  shortest  of  these  survivors  corresponds  to 
^  the  most  likely  discrete  state  sequence. 

The  algorithm  can  be  stated  formally  as  follows: 

INITIALIZATION:  k  *  0 

d[d(0)]  -  dQ 

r[d(03]  =  0 

STEP  k:  d[d(k)  ]  0  <_  d(k)  <_  N 

T[d(k)3  0  <_  d(k)  <_  N 

Compute  r[d(k+l),  d(k)  ]  4  r[d(k)]  [d(k+l)  ,d(k)  ] 

for  all  d(k+l),d(k) 

Find  r[d(k+l)]  -  Min  r[d(k+l)  ,d(k)] 
d(k)6  D 

Store  r[d(k+l)]  and  the  corresponding  d[d(k+l)] 
k  ■  k+1  and  repeat  until  k  =  K. 

STEP  K:  d[d(iQ3  0  <_  d(K)  <_  N 

r[d(K)J  0  <_  d(K)  <_  N 

Find  Min  r[d(K)]  . 

d(K)6D 

/v 

The  corresponding  survivor  d[d(K)]  is  d  .  END. 

The  main  feature  in  this  algorithm  is  that  we  need  only  to  store 

paths  at  each  time  k  ■  0,1,,.., K.  A  straightforward  computation  of 

K  K 

would  have  needed  the  storage  of  N  paths  at  time  K. 

In  practice,  certain  trivial  modifications  are  necessary.  This 


is  the  subject  of  the  following  section. 
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4.3  Actual  Implementation 

When  the  state  sequences  are  very  long  or  infinite,  it  is  necessary 
to  truncate  the  survivors  to  some  manageable  length  Kf.  The  algorithm 
must  decide  on  nodes  up  to  time  k  -  at  time  k.  If  the  truncation 
depth  is  chosen  large  enough,  there  is  a  high  probability  that  all 
time-k  survivors  will  go  through  the  same  nodes  up  to  time  k  -  Kr;  in 
this  case,  the  truncation  costs  nothing.  If  not,  then  we  choose  the 
time  (k  -  node  associated  with  the  shortest  survivor.  If  K?  is 
large  enough,  then  the  effect  on  the  performance  is  negligible. 


In  this  chapter,  we  have  presented  the  Viterbi  algorithm,  which 
is  an  optimal  solution  to  the  problem  of  estimating  the  state  sequence 
of  a  discrete  Markov  process. 

We  next  consider  its  continuous  counterpart,  the  MAP  state  estima¬ 
tion  of  a  continuous  Markov  process. 


4 
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CHAPTER  S 

STATE  ESTIMATION  OF  A  CONTINUOUS  MARKOV 
PROCESS  -  KALMAN- BUCY  EQUATIONS 

In  this  chapter,  we  discuss  the  problem  of  estimating  the  state 
of  a  continuous  Markov  process.  The  Bayesian  formulation  leads  to  a 
recursive  form  for  computing  the  MAP  estimates.  We  will  see  that  the 
concepts  are  identical  to  the  discrete  state.  An  additional  Linear- 
Gaussian  assumption  gives  some  computational  elegance  with  the  Kalman 
Equations . 

S.l  Problem  Statement 

The  system  considered  is  a  continuous  dynamic  system  observed 
during  an  interval  [0,K].  The  system  model  is  the  one  developed  in 
Chapter  2  Section  2.  Here,  at  each  time  k  *  0,1,..., K,  the  discrete 
state  d(k)  can  only  take  one  value,  so  that  the  dependence  in  d(k)  is 
dropped.  Thus,  we  have: 

P(x(k)!xk_1)  •  P(x(k) i x(k>l)}  (5.1) 

P(zOO  |xk,zk-1)  «  P(z(k)|x(k))  (S.2) 

We  have  also  the  following  simplified  equations  for  the  Linear 
Gaussian  assumption: 

x(k+l)  *  F(k)x(k)  ♦  w(k)  (3.3) 


z(k)  »  H(k)x(k)  +  v(k)  , 


(5.4) 
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where  x(0),  w(k)>  v(k)  are  independent  zero-mean  White  Gaussian  Processes 
with: 

E(x(0)xT(0)}  *  PQ 

E{ w(lc)wT(k)}  -  Q(k)  (S.S) 

Ef  vOOvT(k)}  *  R(k)  . 

In  many  estimation  problems,  the  objective  is  to  find,  at  each 
time  k  *  0,1,..., K,  an  MAP  estimate  of  the  instantaneous  state  x(k)  of 
a  system,  based  on  the  available  measurement  data  at  that  time:  this 
is  a  filtering  problem. 

However,  in  some  cases,  a  refinement  of  this  estimation  is  neces¬ 
sary.  For  instance,  in  the  launch  of  a  satellite,  it  is  desirable  to 
determine  the  performance  of  tha  guidance  and  propulsion  systems  during 
the  initial  phases  of  the  flight,  so  that  an  estimation  of  the  state 
traj ectory  has  to  be  done,  using  measurement  data  beyond  the  current 
time.  This  is  a  smoothing  problem,  whose  objective  is  the  MAP  estimation 
of  the  state  sequence  x  *  (x(0) , . . . ,x(K)) ,  based  on  all  available 
measurement  data. 

3.2  The  Filtering  Problem 

As  mentioned  above,  the  objective  in  the  filtering  problem  is  to 
find  for  each  time  k  ®  0,1,..., K,  the  MAP  estimate  of  the  state  x(k) , 
which  is  given  by  x(kjk)  for  which  the  a  posteriori  density  P(x(k)|zk) 
is  maximized.  In  short,  the  problem  is:  for  all  k  *  0,1,... ,K,  find 
$(kfk)  such  that: 

V  x(kj  €  X,  P[*Ck(iOl=k]  1  P[x(k)jzk]  . 

This  estimation  problem  consists  first  in  the  determination  of  the 
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a  posteriori  density  P(x(k)J  zk) ,  then  in  the  maximization  of  this  function 
over  the  state  space  X. 

Using  Bayes'  rule  and  (3.1)  and  (S.2),  we  can  determine  recursively 
the  density  P(x(k)|z  )  from  the  following  relations: 

for  all  k  *  0,1, . . .  ,K  , 

update:  P(x(k)|zk)  »  x  P(x(k)|  zk_1)  (5.6) 

P(z(k)|  z  A) 

normalize:  P(z(k) | zk_1) *  /  P(x(k)izk"1)  x  P(2(k) [ x(k)) dx(k)  (5.7) 

X 

predict:  P(xCk) | zk_1) *  S  P(x(k) j x(k-l))  x  P(x(k-1) | :k'1)dx(k-l)  (5.8) 

X 

initialize:  P(x(0)|z-1)  =  P(x(0))  . 

J- 

Then,  having  computed  P(x(k)jz  ),  we  use  standard  optimization 
techniques  to  determine  the  state  estimate  x(kjk). 

The  main  difficulty  in  this  problem  lies  in  the  explicit  determina- 
tion  of  P(x(k)|z  ),  so  that  it  is  often  impossible  to  determine  the  MAP 

A 

estimate  x(kjk). 

However,  under  the  additional  Linear-Gaussian  assumption  (5.3)- 

k 

(5.S),  the  density  P(x(k)|z  )  is  Gaussian,  and  so  is  completely  defined 
in  terms  of  a  finite  set  of  parameters:  the  conditional  expectation. 
E{x(k)jzk}  and  the  covariance  E((x(k)  -  E(x (k) J zk}) (x(k)  -  E(x(k) )  zk)T)  zk}  . 
This  introduces  a  great  simplification  in  our  estimation  problem,  as  we 
shall  see  in  the  next  section  with  the  Kalman  Filtering  Equations. 
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5.3  The  Kalman  Filtering  Equations 

i  A 

Since  P(x(k)j2  )  is  a  Gaussian  distribution,  the  value  x(k|  k) 

k  k 

for  which  P(x(k)|z  )  is  maximum  is  simply  the  conditional  mean  E£x(k)[z  }. 

In  1960,  Kalman  [2]  derived  a  method  for  computing  recursively 

fc(k[k)  and  the  corresponding  covariance  P(kjk). 

Throughout  this  work,  &(i[j)  denotes  the  optimal  estimate  of  x(i) 

based  on  measurement  data  up  to  time  j,  and  P(i|j)  the  corresponding  co- 

variance  matrix. 

The  Kalman  Filtering  Equations  compute  recursively  £(k[k)  and 


P(k|k)  as  follows: 

£(k|k)  -  x(kik-l)  ♦  K(k)[z(k)  -  H(k)x(k| k-1)]  (5.10) 

fc(k+l[ k)  »  F(k)5(k|k)  (5.11) 

P(k[ k)  *  [I  -  K(k)H(k)3PCk|k-l)  (S.12) 

P(k+l|k)  -  F(k)P(k| k)FT(k)  +  Q(k)  (5.13) 

K(k)  =  P(kjk-l)HT(k)[H(k)P(k|k-l)HT(k)  +  R(k)]"1  (5.14) 

*(0|-1)  «  xQ  PCOj-1)  -  PQ  (5.15) 

F(k),H(k) ,  Q(k) ,  R(k) ,  PQ  are  defined  by  (5. 3) -(5. 5). 


k 

Thus,  the  a  posteriori  density  P(x(k)|z  )  is  a  Gaussian  distribu¬ 
tion,  whose  mean  $(k}k)  and  covariance  P(k|k)  are  simply  determined  by 
(5.10)-(5.1S). 

Another  distribution  of  interest,  as  we  shall  see  later,  is  the 
probability  density  for  the  innovation  process  (:(k)  -  H(k)£(k| k-i) , 
k  »  0,1,..., K}.  This  process  is  a  zeTO-mean  White  Gaussian  Noise  with 

k 

a  covariance  B  equal  to 

Bk  «  H(k)P(k!k-l)HT(k)  ♦  R(k)  (5.16) 
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The  Kalman  Filtering  Equations  provide  directly  the  set  of  estimates 
(2(k}k),  k  *  0,1,..., K}  through  equations  (5. 10) -(5. 15)  by  implicitly  com- 
puting  the  probability  density  P (x Ck) !  2  ).  This  is  the  optimal  recursive 
implementable  solution  to  our  filtering  problem. 

We  next  discuss  the  smoothing  problem  which  uses  additional  measure¬ 
ment  data  to  compute  the  set  of  estimates. 

5.4  The  Smoothing  Ptoblem 

The  objective  in  the  smoothing  problem  is  to  find  an  MAP  estimate 

K 

of  the  state  sequence  x  ■  (x(0) , . . .  ,x00)  based  on  all  available  measure- 

Y 

ment  data  z  »  (z (0) , • . • , z (K) ) . 

kI  k 

The  solution  to  this  problem  is  given  by  the  sequence  X  1  * 

{ £(k|  K) ,  k  *  0,1,..., K}  for  which  PCx^z*)  is  maximum.  In  short,  the 
problem  is: 

find  $K!K  such  that  V  xK  6  XK,  PC**!*!**)  >  PCx^z*)  . 

k!  k 

A  straightforward  computation  of  £  1  would  consist  first  in  the 

K  K  K  1C 

determination  of  P(x  [z  ),  then  in  the  optimization  of  P(x  ]:  ]  over 

the  super-space  X  . 

In  fact,  because  of  the  raemoryless  and  Markov  properties  (5.1)- 
(S.2),  we  shall  see  that  this  estimation  problem  is  equivalent  to  a 
"continuous"  Viterbi  algorithm,  which  requires  much  less  computation. 

For  this  purpose,  let  us  define  the  following  function  at  time  k: 

P(x(k)|zk)«  Sup  P(xk|zk)  (5.17) 

xk"lSXk"1 

We  denote  by  X[x(k)]  the  corresponding  sequence  terminating  at 
x(k) .  It  is  the  equivalent  of  the  survivor  d[d(k)]  in  the  Viterbi  algo- 
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rithm,  in  the  sense  that  the  most  likely  sequence  k  1  must  begin  with 
the  sequence  kikCkj  K)  3  all  'K  *  0,1,..., K.  This  is  true  because  of 
the  following: 

V  xk,  P(xkj  2k)  <.  P(x(k)j zk)  *  Pfk(x(k))|  2k) 

«>  V  xk  €  XK,  V  2(k+l),...,z(K), 

P(xKf  2KD  1  PC^CxCk)).  x(k+l)  , .  •  •  ,x(K)|  zk)  using  (S.l)-(5.2) 

->  Sup  P(xK|zK)<  Sup  P[k[x(k)  ]  ,x(k*l) ....  ,x(K)|  zK] 

x^GX*  "  x(k),...,x(K)GX 

=»>  Sup  PCxV)  *  Sup  P[£(x(k)) ,x(k+l), . . . ,x(k)|  z^] 

x*€XK  x(k),...,x(K)GX 

(5.18) 

K  K 

But  P(x  |  z  )  is  maximum  for{k(k|K),  k  *  0,l,...,k}.  Assuming  a 

K  K 

single  maximum  for  the  density  P(x  |s  ),  we  have:  k[£(kjK)]  »  (k(Q|k),..., 
£(k[  K) )  using  (4.18).  Q.E.D. 

Thus,  fox  all  k  ■  0,1,...,K, 


*[X(k|fC)]  «  (X(0|K),...,X(k|K)) 


(5.19) 


k 

that  is,  we  need  only  to  store  the  sequence  £[x(k)]  and  P(x(k)|z  )  at 
time  k. 

Using  (5.1)-(5.2),  we  can  determine  P(x(k)jz  )  recursively  from 
the  following  equation: 


p(»P0|ik)  .  ,  Sup  P(x(k)[x(k-1))  x  P(x(k-l)l2k_1) 

P (2 (lc)  [  2*”  j  x(k-l)€X 


(5.20) 

where  P(2(k)|zk_1)  is  defined  by  (5.7)-(5.8)  and  P(x(0)iz°)  •  P(x(0)j  :(0)) 


is  the  initial  condition. 


4: 


From  this,  we  c*n  compute  the  functions  i'x.k  ■  recursively.  The 

recursion  proceeds  as  follows:  3t[xvk>  ]  is  extended  oy  zr.t  tiae  urit  and 

the  corresponding  P(x(k) ,x(krl);  zk**)  *  - » -  >-■* 4 1 — — *  P(x(k*L) , x(k))  * 
1c  PC=(k*l),  :  , 

P(x(k)[z  )  is  computed  for  all  x(k),  x(k*l),  using  the  received  aeasure- 

k  ^1 

ment  z(k+l)  at  tiae  k*l.  Then  the  maximization  of  A(x(krl) ,x(k)|  z  ) 
over  x(k)  €  X  leads  to  P(x(k>l}|  z  )  according  to  (5.20).  The  corres- 
ponding  ^[xCk+1)  ] ,  that  maximizes  P(x(k*l)  ,x(k)j  z  ),  is  used  for  the 
computation  of  X(x(k+1)],  since  we  have: 


5t[x(k+l) ]  -  fltp^Cxfkn))],  x<VU)  (5.21) 

At  time  K,  the  maximization  of  P(x(K)!zK)  over  X  yields  5t(K|  K) 

jf!  Y  ]( 

and  subsequently  the  MAP  estimate  X  1  of  the  sequence  x  ,  since,  accor¬ 
ding  to  (5.19),  we  have: 


S[*(K|K)J  -  (XCk|K).  k  -  0,1,..., K}  . 

The  algorithm  can  be  stated  formally  as  follows: 

STEP  0:  P(x(0)|z°)  -  PCx(0)|z(0)) 

SfjcCO)]  -  xCO) 

STEP  k  P(x(k)[zk)  ,  x(k)  6  X 

fc(x(k) ]  ,  x(k)  €  X 

Compute  P(x(k),x(k+l)]zk+1)  -  x 

P(z(k*l)izK) 

P(x(k+1) i x(k))  x  P(x(k)|zk)  , 
for  all  x(k) ,  x(k*l) . 

Find  £(x(k+l) | zk+1)  ■  Sup  P(x(k) ,x(k*l) j zk+1) 

x(k)6X 

•  k*l 

Store  P(x(k-*-l)jz  )  and  the  corresponding  X^(x(k+1)) 
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Then  x[x(k+l) ]  *  (afl^CxCk+l))],  x(k+l)) 
k  *  k+1  and  repeat  until  k  =  K. 

STEP  K:  PCxOOlA  ,  x(K)  6  X 

&[x(K)]  ,  x(lQ  S  X 

Find  Sup  P(x(K)jzK).  The  corresponding  x(K)  is  equal  to 
x(K)6X 

$(K|K)  and  { £(kj  K) ,  k*0,l,...,K}  **[*(K(K)]  . 

END. 


This  algorithm  can  be  considered  as  the  conceptual  extension  of 
the  Viterbi  algorithm  to  the  continuous  case:  at  each  time  k,  we  need 
only  to  store  the  "survivor"  X[x(k)]  and  the  function  P(x(k)',z  ). 

Of  course,  the  main  difficulty  of  this  method  lies  in  the  explicit 
determination  of  P(x(k)|z  )  and  S[x(!t)].  In  the  general  case,  such  a 
determination  is  usually  impossible.  However,  as  in  the  filtering  prob¬ 
lem,  the  Linear  Gaussian  assumption  introduces  a  great  simplification  in 
our  smoothing  problem  with  the  Linear-Gaussian  Smoothing  Equations. 

5.5  The  Linear-Gaussian  Smoothing  Equations 

Under  the  Linear-Gaussian  assumption  (5. 3) -(5. 5),  the  a  posteriori 
probability  distribution  P(xK|  z*~)  is  Gaussian.  Thus,  the  MAP  estimate 
of  xK,  which  is  the  state  sequence  K  ■  {  X(k|  K) ,  k  «  0,1,..., K}  for 

VS  • 

(axpoiz"},  k  -  0 . ,k). 

Meditch  in  [3]  gives  a  method  for  computing  recursively  the  set 
of  smoothed  estimates  (£(k|K),  k  ■  0,1,..., X}  and  the  corresponding  co- 
variances  {P(kjK),  k  *  0,1,..., K}  . 


which  P(xK|  z*S  is  maximum,  is  simply  the  conditional  mean  E{  x 
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The  Linear-Gaussian  Smoothing  Equations  are  the  following:  for 
k  »  K-l,K-2, . . . ,0, 

x(k|K)  *  *(k[k)  +  L(k)  [$(k+l|  K)  -  £(k+l|k)]  (5.22) 

k(K|  K)  is  the  boundary  condition  for  k  »  K-l 

L(k)  *  P(k[ k)FT(k)P-1(k+l|  k)  (5.23) 

P(k|K)  »  P(k|  k)  +  L(k)  [P(k+l[  K)  -  P(k+l[ k) ]LT(k)  (5.24) 

X(k|k),  XCk+l[k),  P(k|k),  P(k+l|k)  are  given  by  the  Kalman  Fil¬ 
tering  Equations  (5.10)-C5.1S) . 

The  filtering  equations,  as  we  saw  in  Section  3,  compute  recur¬ 
sively  forwards_Jji__yjne  the  set  of  estimates  (£(kjk),  k  ■  0,1,..., K}  and 

4 

the  corresponding  covariances  {P(kjk),  k  *  0,1,...,K}.  Then  the  smoothing 

equations  (5. 22) -(5. 24)  compute  recursively  backwards  in  time  the  set  of 

smoothed  estimates  {*(k[K),  k  ■  K-l.K-2, . . . ,0}  and  the  corresponding 

covariances  (P(k|K),  k  *  K-1,...,0}. 

Thus,  in  the  Linear  Gaussian  case,  the  smoothing  equations  simply 

refine  the  set  of  estimates  given  by  the  filtering  equations. 

The  set  of  smoothed  estimates  is  computed  directly  without  de- 

K 

termining  explicitly  the  functions  TC-xCk)  [  z  )  for  k  *  0,1,..., K.  Later 

k 

in  this  work,  an  expression  for  7(x(k)|z  )  will  be  needed.  The  probabi¬ 
lity  distribution  P(xk|zk)  is  equal  to  N[xk  -  xk|k,pk|kj.  It  can  j)e 
shown  that  P(x(k)|zk)  ■  Sup  P(xklzk)  is  equal  to: 

P(x(k)  |  zk)  N[x(k)-*(k|k),  P (k |  k)  J  ,  (5.25) 

^C2n)nk|pk  k[ 

where  li(kjk),  P(k[k)  are  given  by  the  Kalman  Filtering  Equations  (5.10)- 

k  i  k  k 

(5.15)  and  P  1  is  the  covariance  of  the  overall  state  sequence  x  .  (See 


4S 


Fl&URE  S.A 


x(A|t)  iffcfK)  3 

©  P(z(t)l**),  N[x(l).£((Llk),  p(4 JK)  ] 

©  P(*f4)|z*j»  N[x(*)-£(*H),  P(i /A)] 


©  PfzWlz4)g\|(^1PM,  P(,(*)|»*J 


iu)  P' 
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Fig.  5.1). 

5.6  Summary 

In  this  chapter,  we  have  discussed  the  MAP  continuous  state  esti¬ 
mation  problem.  The  difficulty  in  the  filtering  problem  lies  in  the 
determination  of  the  density  P(x(k)|z^).  However,  under  a  linear-Gaussian 
assumption,  the  set  of  estimates  (x(k|k),  k  ■  0,1,..., K)  can  be  computed 
directly. 

The  difficulty  in  the  smoothing  problem,  which  can  be  viewed  as 
a  "continuous"  Viterbi  algorithm,  lies  also  in  the  determination  of  some 
functions  PCx(k)jz^)  and  x[x(k)].  Again,  the  Linear-Gaussian  assumption 
introduces  a  great  simplification,  in  the  sense  that  a  direct  recursive 
computation  of  the  set  of  smoothed  estimates  (x(k[K),  k  *  K-1,...,0)  be¬ 
comes  possible.  Besides,  this  set  of  estimates  is  simply  obtained  by 
"refining"  the  set  of  filtered  estimates. 

Now,  having  presented  the  Viterbi  algorithm  in  Chapter  4  and  the 
MAP  continuous  state  estimation  problem  in  this  chapter,  we  discuss  in 
the  next  chapter  the  combined  MAP  hybrid  state  estimation  problem. 
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CHAPTER  6 

STATE  ESTIMATION  OF  A  HYBRID  MARKOV  PROCESS 


In  this  chapter,  we  consider  a  single  hybrid  system.  The  purpose  of 
this  part  is  to  obtain  an  algorithm  for  estimating  the  hybrid  state  sequence, 
that  keeps  the  amount  of  computation  at  a  reasonable  level.  The  combined 
MAP  hybrid  state  estimation  algorithm  in  a  blend  of  the  Viterbi  algorithm 
of  Chapter  4  and  the  MAP  continuous  state  estimation  algorithm  of  Chapter 
S.  Depending  on  two  possible  choices  of  objective  for  this  problem, 
there  are  two  different  formulations,  as  we  shall  next  see, 

6.1  Problem  Formulations 

The  model  for  the  hybrid  system  considered  is  the  one  of  Chapter 
2  Section  2.  The  system  is  observed  during  an  interval  [0,K].  The 

K 

problem  is  to  estimate,  from  the  available  measurement  data  z  =  C-C°) » • • • . 

z(K)),  the  hybrid  state  sequence  s  *  Cs(0)  , . .  •  ,s(K)) . 

K  K 

The  MAP  estimation  of  s  based  on  z  is  one  possible  objective. 

K  jf 

In  this  case,  the  problem  is  to  find  s  =  id ,x j  such  that: 

K  _  VK  „  .K  _  _K  _,0K  JK,  K.  .  K  ,K,  K. 

V  x  ex,  V  d  e  D  ,  P($  ,d  1  >  P(x  ,d  !  z  )  . 

However,  in  some  cases,  more  importance  is  given  to  the  estimation 
of  the  discrete  state  sequence  d  .  For  instance,  in  failure  detection 
problems,  when  the  system  has  sufficient  hardware  back-up  capabilities, 
we  could  be  interested  only  in  detecting  the  failures,  without  directly 
desiring  to  estimate  the  continuous  state  of  the  system.  Similarly,  in 
some  surveillance  problems,  we  can  be  more  interested  in  the  number  of 


targets  in  a  certain  area,  than  in  their  precise  locations. 


In  these  cases,  we  first  perform  a  MAP  estimation  of  the  discrete 

K  K 

state  sequence  d  :  so,  the  problem  is  to  find  d  such  that: 

V  dK  6  DK  ,  P(dK|  2*)  >  P(dK|  2*) 

Then,  if  we  are  still  interested  in  the  continuous  state  of  the 

U 

system,  an  MAP  estimation  of  the  sequence  x  can  be  done,  conditioned 

AjJ 

on  d  ;  this  problem  is  simply  the  continuous  state  estimation  problem 
discussed  in  Chapter  S.  In  this  objective,  we  minimiae  the  probability 
of  an  error  in  estimating  the  discrete  state  sequence  d  . 

Both  objectives  do  not  yield  necessarily  the  same  estimates 

The  first  approach  leads  to  a  better  estimate  of  the  pair  (dK 
and  the  second  one  to  a  better  estimate  of  the  discrete  sequence  dK. 

In  both  cases,  a  graph  can  be  associated  with  our  Markov  process, 
as  in  the  Viterbi  algorithm.  In  this  graph,  each  node  corresponds  to  a 
distinct  discrete  state  d(k)  of  the  system,  and  each  branch  represents 
a  transition  to  some  new  discrete  state  at  the  next  instant  of  time. 

The  graph  begins  at  a  known  discrete  state  dg.  Thus,  to  every  possible 
discrete  state  sequence  d  ,  there  corresponds  a  unique  path  through 
the  graph,  and  vice  versa. 

In  the  next  two  sections,  we  consider  both  approaches,  and  derive 
for  each  case  some  pruning  rules  in  the  graph  that  are  available  to  re¬ 
duce  the  amount  of  computation. 

6.2  Discrete  State  Estimation  Algorithm 

As  mentioned  above,  under  this  objective  we  have  to  perform  a  MAP 

K  K 

estimation  of  the  discrete  state  sequence  d  :  find  d  such  .that: 
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V  dK  €  DK  ,  P(dK|zK)  >  P(dK[zK) 

£ 

A  straight forward  computation  of  d  is  combinatorial:  find  the 

K  K 

most  likely  state  sequence  d  among  N  state  sequences.  A  method  has 
to  be  found  for  decreasing  the  amount  of  computation. 

This  problem  is  formally  equivalent  to  the  problem  of  finding 

It 

the  shortest  path  in  the  above  graph,  if  we  define  the  length  of  X(d  ) 

jt 

of  a  path  corresponding  to  the  state  sequence  d  by: 

X(dk)  »  -in  P(dk,zk)  (6.1) 

However,  the  solution  to  this  shortest  path  problem  is  not  as 

simple  as  for  the  discrete  state  case.  The  Viterbi  algorithm  uses  the 

Markov  and  memoryless  properties  to  keep  only  the  N  survivors  u[d(k) ] 

at  each  time  k  *  0,1,..., K.  But  here,  the  memoryless  property  does  not 

i  k- 1  k 

hold  anymore:  the  probability  distribution  P(z(k)|z  ,d  )  depends  on 

v 

the  whole  discrete  state  sequence  d  and  on  all  the  past  measurements 
k-1 

z  through  the  continuous  states  of  the  hybrid  system.  Thus,  the 
k  k 

condition  X(d^)  ^X(d2),  where  d^Ck)  ■  d2(k),  is  not  sufficient  anymore 

k  k 

for  discarding  d2  at  time  k.  As  later  measurements  are  received,  d2 

'  K 

can  well  be  a  part  of  the  most  likely  sequence  d  . 

A  stronger  condition  is  needed  for  discarding  a  path  corresponding 
v 

to  d  at  time  k,  without  increasing  the  probability  of  an  error  in  esti- 
mating  the  whole  sequence  d  .  Th^s  condition  is  a  consequence  of  the 
following  theorem: 

k  k 

Theorem  1:  Let  dj  and  d2  be  two  discrete  state  sequences  such  that 
dx(k)  -  d2(k).  If  for  all  x(k)  S  X,  P(x(k),  dk[:k)  <_  P(x(k),  dkj:k), 
then  for  all  z(k+l) , . . . ,z(K)  and  for  all  d(k+l) , . . . ,d(K) ; 


4 
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P(dk!  zK)  <_  PCdJls'S  where  d*  -  (dk,d(k+l) , . . .  ,d(K))  and  d*  =  (dk,d(k+l) , . . . 
dffl). 

Proof:  The  proof  of  this  theorem  is  straightforward  using  the  following 

ki  k 

recursive  relation  on  P(x(k),d  |z  ): 

P(x(k),dk|  zk)  »  ?.C2.Q?)J^k^,xCk)j  x  pCxCk)  >dk]  2k‘l)  (6.2) 

P(z(k) 1 zK  j 

PCr(k)|zk"1  »  Z  /  PCxOO.dV'’1)  ^ 

{dk}  X 

P(z(k)  |x(k) ,  d(k))dxCk)  (6.3) 

P(x(k),dk[zk"1)  »  /P(x(k),d(k)|x(k-l),dCk-l)  x 
X 

P(x(k-1) ,dk”1 | zk"i)dx(k-l)  (6.4) 

0  -1  A 

P(x(0),d  |z  )  *  P(x(0))  is  the  initial  condition  (6.5) 

Thus,  simply  using  (6.2),  if  V  x(k),  P(x(k) ,dk| zk)  <_  P(x(k)  ,dk  jzk) 
where  dj(k)  *  d^(k),  then  V  x(k+l),  V  d(k+l) ,  V  z(k+l), 

P(x(k+l),dk+1|zk+1)  <_  P(x(k+1) ,dk+1 jzk+1) .  And  proceeding  by  induction: 

V  x(K),  V  d(k+l)  ....  ,d(K)  ,  V  z(k+l)  , . . .  ,z(K)  ,  P(x(K),  d*|z*)  < 

P(x(k)  ,di^ jz*) .  Finally,  by  integrating  on  x(K) :  P(d^  1  z*S  P Cd^  i  • 

Q.E.D. 

This  theorem  means  that  we  can  compare  two  discrete  state  se¬ 
quences  at  time  k,  terminating  at  the  same  node  d(k) :  if  at  this  time, 
all  the  continuous  states  are  less  likely  on  one  sequence  and,  if  after 
that  time  we  follow  a  same  common  path,  then  this  trajectory  will  remain 
less  likely,  whatever  the  future  measurements  are.  Therefore,  we  can 


eliminate  it  immediately  from  consideration.  Thus  we  have  the  first 
pruning  rule  in  our  graph: 


k  k 

Pruning  Rule  R1 :  Let  d^  and  d^  be  two  discrete  state  sequences  such  that 
dx00  -  d2Ck).  If,  for  all  x(k) ,  P(x(k),  dk!zk)  <  P(x(k) ,dk| zk) .  (6.6), 
then  the  whole  sequence  d^  can  be  immediately  discarded  at  time  k,  with- 
out  increasing  the  probability  of  an  error  in  estimating  d  . 

The  implementation  of  this  pruning  rule  is  not  as  simple  as  in 
the  Viterbi  algorithm.  Besides,  it  is  rarely  possible  to  order  each 
discrete  state  sequence  using  (6.6):  most  of  the  time,  more  than  one 
'’survivor"  terminating  at  a  given  d(k)  remain  after  having  used  this 
pruning  rule.  Therefore,  we  would  like  to  extend  the  possibilities  of 
eliminating  some  sequences.  This  is  the  subject  of  the  following  the¬ 
orem,  which  allows  comparisons  between  one  sequence  and  a  set  of  se¬ 
quences  : _ 

k  k  k 

Theorem  1 ' :  Let  dg,d^,...,dr  be  (r+1)  discrete  state  sequences  such 

that  d_(k)  *  d.  (k)  *  ...  d  (k) .  Let  (a.),,.-,  be  r  real  positive 
u  i  ^  r  i  iv  is  r 

Y  —  k  k 

coefficients  such  that  E  »  1.  If  for  all  x(k)  G  X,  P(x(k),  dQ|z  )  <_ 
r  .i»l 

E  a  P(x(k) ,d. (k)  { z  ) ,  then  for  all  r(k+l) , . . . ,z(K) ,  for  all  d(k+l),... 

i*l  *  ^ 

d(K) ,  PfdlJlz’S  <  Sup  P(df|zK),  where  for  all  i  »  0,1,..., r,  df  * 
v  l<i<r 

(df,d(k*l),...,d(k)>“ 

Proof:  The  proof  of  theorem  1'  is  very  similar  to  the  one  of  theorem 
1  and  is  omitted. 


We  notice  that  the  set  of  a. 's  is  arbitrary  as  long  as  a.  >  0 
r  1  1  “ 

and  E  dy  ■!,.  A  procedure  for  choosing  (<x)  could  be: 
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Max  Inf  Z  a  P(x(k)  ,dk|  zk)  -  P(x(k)  ,dkj  zk) 
a  x(k)GX  i=l 

r 

s.t.  a.  >  0  and  I  a.  a  i 
1  —  .  .  i 

i»l 

If  this  Max/Min  is  positive,  then  we  have:  P(d^| z*S  <  Sup  P(d!f|z 

Anyway,  from  Theorem  1*  we  deduce  immediately  the  following  pru¬ 
ning  rule: 

V  V  V 

Generalized  Pruning  Rule:  GR1:  Let  dg,d^,...,dr  be  (r+1)  discrete  state 

sequences  such  that  dg(k)  »  ...  »  dy(k) .  If  it  exists  r  positive  coef- 

r 

ficients  ct.  such  that  Z  a.  =  1  and 
1  i»l  l 

V  x(k)  G  X  ,  P(x(k) ,  dkjzk)  <  Z  a.P(x(k),dJlzk)  ,  (6.7) 

i=l 

then  d*  can  be  immediately  discarded  at  time  k,  without  increasing  the 

K 

probability  of  an  error  in  estimating  d  . 

GR1  is  a  more  performant  pruning  rule  than  Rl,  in  the  sense  that 
more  state  sequences  can  be  discarded  at  a  given  time  using  that  rule 
(it  reduces  to  Rl  when  r  =  1) ;  but  it  is  also  more  difficult  to  implement. 

Using  Rl  or  GR1,  we  are  now  capable  of  discarding  some  sequences 
during  the  formation  of  the  graph.  Of  course,  since  these  rules  are 
based  on  functional  inequalities,  more  than  one  "survivor"  will  often 
remain  for  every  sequence  terminating  at  a  given  node.  At  time  K,  we 
choose  d  among  the  remaining  state  sequences.  This  is  a  substantial 
progress  on  the  straightforward  computation. 

The  algorithm  can  be  stated  formally  as  follows:  letS  be  the 
set  of  remaining  sequences  at  time  k. 


ALGORITHM  1 


STEP  0:  J  2?  =  dQ 

I  P(x(0) ,  d°I2°)  *  PCx(0)iz(0),d0) 

STEP  k:  f^k-l 

(  P(x(k-1)  Jdk“1tzk_1)  for  all  dk_1  G^"1,  x(k-l)  6  X 


Compute  P(x(k) ,dk| zk)  using  (6. 2) -(6. 5),  for  all  x(k)  G  X, 

k  k-1 
dKGSI  1  x  D. 

Apply  GR1  to  dk:  for  all  dk  GS^"1  *  D  -  (dk), 

YES 

Max  Inf  Za.  P(x(k) ,dj| zk)  -  P(x(k) ,dk| zk)  >  0 

a.  x(k)GX  i  1 

s.t.  ai  >_  0  and  la^  *  1 

YES:  Delete  dk 

NO:  Store  dk  in  2^ 

k  *  k+1  and  repeat  until  k  ■  K 


STEP  K: 


(  P(x(K) ,  dK|zIC)  for  all  x(K)  6  X,  dkSSK 


Compute  P(dK|zK)  -  /  P(x(K),  dK|  z’SdxOO  ,  for  all  dK  6  2^. 

X 


Find  Sup  P(d|z 
dK€2>K 


jr  K 

) .  The  corresponding  d  is  equal  to  d 


END. 


If  we  are  also  interested  in  an  estimate  of  the  continuous  se- 
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v 

quence  x  ,  we  have  to  combine  Algorithm  1  with  the  MAP  continuous  state 
estimation  algorithm  or  Chapter  5,  and  keep  also  track  at  time  k  of 
P(x(k),dkizk)  and  x[x(k),dk]. 

Of  course,  the  difficulty  in  algorithm  1  lies  in  the  explicit  de- 

V  i  V 

termination  of  P(x(k),d  \z  )  and  in  the  application  of  GR1.  The  Linear- 
Gaussian  assumption  introduces  a  great  simplification  in  these  computa¬ 
tions,  as  we  shall  see  in  Chapter  7.  We  next  discuss  the  second  objec¬ 
tive  of  a  hybrid  state  estimation  problem. 

6.3  Hybrid  State  Estimation  Algorithm 

The  objective  in  this  part  is  to  find  the  MAP  estimate  of  the 

K 

hybrid  state  sequence  s  based  on  all  available  measurement  data.  The 

K  ''K  KI  K 

solution  to  this  problem  is  given  by  §  *  (d  ,£  1  )  such  that: 

V  dK  €  DK-,  V  xK  6  XK  ,  P(dK,xK^  K[  z*S  >  PCdK,xK|zK)  . 

Y 

This  straightforward  computation  of  s  is  performed  as  follows:  first, 
on  each  sequence  d  ,  the  problem  is  to  find  the  continuous  state  se¬ 
quence  x^Cd1^  =  (xQ|  j^Cd^  , . .  .,xKj  jjCd1^)  for  which  PCxK]dK,zK)  is 

maximum:  this  is  the  smoothing  problem,  conditioned  on  d  ,  which  has 

K 

been  presented  in  Chapter  5.  Then,  the  estimation  of  d  consists  in 

''K 

finding  the  sequence  d  such  that: 

V  dK  G  DK,  P(dK,  JcKI  K(31C)  |  z*S  >  P(dK,jc^  K(dK)  |  ,*)  (6.8) 

''K 

The  evaluation  of  d  is  strictly  combinatorial:  find  the  sequence  d 
verifying  (6.8)  among  N  state  sequences.  So,  as  for  the  Viterbi 


algorithm,  we  have  to  find  a  way  for  reducing  the  amount  of  computation. 
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This  problem  is  equivalent  to  the  one  of  finding  the  shortest  path 

jc 

in  the  graph  considered  in  Section  1,  if  we  define  the  length  X(d  )  of 

1c 

a  path  corresponding  to  the  state  sequence  d  by: 

X(dk)  ■  -  In  Sup  P(xk,dk|zk)  (6.9) 

xk6Xk 

lc  1c 

As  for  the  first  objective,  the  relation  X(dj)  <_X(d2)  where  d^(k)  »  d2(k), 

k 

is  not  sufficient  for  discarding  d2  at  time  k,  since  it  does  not  imply 
that  for  all  z(k+l) , . . . ,z(K) ,  and  for  all  d(k ■*■!),...  ,d(K) ,  A Cd^j  z*S  i 

X(d*|zV 

A  stronger  condition  is  needed  for  discarding  optimally  a  path 

corresponding  to  a  state  sequence  d  at  time  k.  This  condition  is  based 

•  ki  k  k  k i  k 

on  the  recursively  computable  functions  P(x(k),d  |z  )  =  Sup  P(x  ,d  |z  ) 

xk‘1exk_1 

It  is  a  consequence  of  the  following  theorem: 
k  k 

Theorem  2:  Let  dj  and  d 2  be  two  discrete  state  sequences  such  that 

dx(k)  »  d2(k).  If,  V  x(k)  €  X,  P(x(k) ,dk| zk)  <.  P(x(k) ,dk| zk) ,  then 

V  z(k+l),...,z(K),  V  d(k+l),...,d(K),  Sup  P(xK,df[  z*)  <_ 

xkgXK 

Sup  PCx^d^lz^  where  dk  *  (dk,d(k+l) , . . .  ,d(K))  and  dk  * 
xK6XK 

(dk,d(k+l),...,d(K)). 

Proof:  The  proof  of  this  theorem  is  straightforward  using  the  following 

k  k 

recursive  relations  on  P(x(k),  d  |z  ), 

p(x(k)  ,dki  zk) .  x 

P(=W!z  1 

Sup  P(x(k)  ,d(k)  j x(k-l)  ,d(k-l) )  x  P(xCk-l)/'1 1  z^1)  (6.10) 
x(k-l)GX 
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PCsClOU14'1)  is  defined  by  (6.3)-(6.4). 

P(x(0) ,d^l z^)  »  P(x(0)jz(0),  d^)  is  the  initial 

condition  (6.11) 

Now,  simply  using  (6.10),  we  have: 

V  x(k)  G  X,  P(x(k),dk|zk)  <_  P(x(k),dk!zk)  where  d^k)  »  d2(k) 

*>  V  x(k+l) ,  V  z(k •*•]),  V  d(k+l),  P(x(k+l)  ,dk+1lzk+1)  <_ 

P(x(k+l),dk+1| zk+1) 

Proceeding  by  induction,  we  have: 

Vx(K) ;  Vz(k+l),...,z(K);  V  d(k+l) , . . . ,d(K) , 

P(x(K),dk|  z*)  <  P(x(K),d![|zK) 

->  V  z(k+l),...,z(K);  V  d(k>l),...,d(K)  , 


As  fOT  the  first  objective,  we  can  compare  two  sequences  termi- 

-  k  k 

nating  at  the  same  state  dCk) .  The  comparison  is  based  on  P(x(k),d  |z  ) 

k  lc  K 

instead  of  P(x(k),d  |z  j,  since  to  delete  an  hypothesis  d*  at  time  K,  we 


P(dkl z*). 


From  Theorem  2,  we  deduce  immediately  the  first  pruning  rule  in 
our  graph. 
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k  k 

Pruning  Rule  R2;  Let  d1  and  d^  be  two  state  sequences  such  that  d^(k)  * 

d2(k).  If,  for  all  x(k)  €  X,  P(x(k) ,dk| zk)  <_  P(x(k) ,d£|  2k) ,  (6.12) 

]{ 

then  d^  can  be  immediately  discarded  at  time  k. 

Since  it  is  not  possible  to  order,  using  (6.12),  all  the  sequences 
terminating  at  a  given  node  d(k) ,  most  of  the  time,  more  than  one  •'survi¬ 
vor"  terminating  at  d(k)  will  remain.  The  possibilities  for  deleting  a 
sequence  are  increased,  if  we  compare  this  discrete  sequence  with  a  set 
of  discrete  sequences  terminating  at  the  same  node.  This  is  the  subject 
of  the  following  theorem: 

k  k  k 

Theorem  2 1 :  Let  d^,d^, . . . .d^  be  (r*l)  discrete  state  sequences  such  that 

dQ(k)  -  ...»  dr(k).  If,  for  all  x(k)  6  X,  P(x(k) ,dkj zk  <_ 

Sup  P(x(k),dkjzk),  then  Vd(k*l) , . . . ,d(K) ,  V  z(k+l) , . ,z(K) , 
l<i<r 

Sup  ?ix*.d*\z*)  <_  Sup  Sup  PCxViz*)  where,  for  all  i  *  0,l,...,r, 
xK6XK 

d£  »  Cdk,d(k+l),...,d(K))  . 

Proof:  The  proof  is  similar  to  the  one  of  Theorem  2  and  is  omitted. 


The  generalised  pruning  rule  2  is  a  direct  consequence  of  this 
theorem: 


Generalized  Pruning  Rule  GR2:  Let  d^,.,.,dr  be  (r+1)  sequences  such 

that  dQ(k)  »  ...  »  dy(k).  If  for  all  x(k)  G  X,  P(x(k)  ,dj|  z5*)  <_ 
k  k  k 

Sup  P(x(k) ,d?| z  ) ,  then  d*  can  be  immediately  discarded  at  time  k. 

K i<r  1  u 

Using  R2  or  GR2,  we  can  discard  some  sequences  during  the  forma- 

jf 

tion  of  the  graph.  At  time  K,  we  choose  d  among  the  sequences  by: 


(6.13) 
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for  all  the  remaining  sequence  d  .  The  estimate  of  the  hybrid  state  se- 

.  -K~K(K  "1C,  -K[K*K,  ..  ,  .  .  u.n  .. 

quence  is  given  by  s  *  (d  ,x  j).  x  (d  ) ,  which  is  the  MAP  esti- 

K 

mate  of  x  conditioned  on  d  ,  is  determined  by  the  continuous  state 

estimation  algorithm  of  Chapter  5.  So,  by  keeping  track,  at  each  time 

_  k  I  lc 

k,  not  only  of  the  functions  PCxQO,d  |z  5,  but  also  of  the  survivors 
a  k  k 

x[x(k),d  ]  for  all  the  remaining  sequence  d  ,  we  will  be  able  to  perform 
at  time  K  the  complete  hybrid  estimation. 

The  algorithm  can  be  stated  formally  as  follows:  let  2^  be  the 
remaining  sequences  at  time  k. 


ALGORITHM  2 

STEP  0:  fgP  =*  dQ 

\  PCxC0},d°|z°)  »  P(xC0)  |  z(0)  ,dQ) 

V  j?lx(0),d°]  »  xCO) 

STEP  k:  (&~l 

\  PCxCk-lbd*'1!  zk_1)  for  all  x(k)€X,  dk_1  G^"1 

v  £[x(k-l)  .d15*1] 


Compute  P(x(k),dk|  zk)  and  $[x(k),dk]  for  all  x(k)  6  X, 
dk  6  S^"1  x  D,  using  C6.10) . 


Apply  GR2  to  dk:  for  all  dk  G^*1  *  D  -  {dk}, 

V  x(k)  6  X,  ?Cx(k),dk|zVfS  Sup  P(x(k) ,dk|  zk) 

No  i 

lr 

YES:  Delete  d 
NO:  Store  dk  in2>k 

k  ■  k-t-1  and  repeat  until  k  *  K 
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STEP  K: 

J  P(.xCK)  ,dK|zK)  for  all  xOO  6  X,  dK  6  Q$L 
(*[x(iq,dK3 

Compute  pc^VS,^/)  .  Sup  P(x(K),dK]  z*) 

x(K)6X 

Find  Sup  P(£^  ^(d*S  ,d^|  z^)  .  The  corresponding 

dKesK 

(3K,XK!  !<(dK))  is  equal  to  kK.  END. 


The  difficulty  in  an  actual  implementation  of  Algorithm  2  lies  in 
an  explicit  determination  of  P(x(k) ,dkj z^)  and  x[x(k),dk],  in  the  storage 
of  these  functions,  and  in  the  application  of  GR2,  which  is  based  on  a 
functional  inequality.  Once  again,  the  Linear-Gaussian  assumption  intro¬ 
duces  a  great  simplification  as  we  shall  see  in  Chapter  7. 

6.4  Summary  and  Conclusion 

In  this  chapter,  we  have  discussed  the  state  estimation  of  a 
hybrid  Markov  process.  Depending  on  the  application,  we  can  be  interested 
either  in  estimating  the  hybrid  state  sequence,  or  in  estimating  only  the 
discrete  state  sequence  of  this  system.  In  both  cases,  we  have  shown  the 
formal  equivalence  between  our  estimation  problem  and  a  shortest  path 
problem  in  a  certain  graph.  Then  we  have  derived  some  pruning  rules, 
based  on  the  distributions  P(xtk),dk|  2k)  or  P(x(k),dk|  zk) ,  that  are 
used  for  the  deletion  of  some  discrete  state  sequences  during  the 
formation  of  the  graph.  This  adaptive  estimation  algorithm  should 
decrease  substantially  the  amount  of  computation  required  at  time  K. 


In  addition  to  an  explicit  determination  of  the  above  distributions, 
there  is  also  the  problem  of  implementing  in  a  simple  way  these  pruning 
rules.  This  can  be  done  under  an  additional  Linear-Gaussian  assumption, 
as  we  shall  see  in  the  next  chapter. 
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CHAPTER  7 

STATE  ESTIMATION  OF  A  HYBRID  GAUSS-MARKOV  PROCESS 


The  purpose  of  this  chapter  is  to  apply  the  results  of  Chapter  6  to 
the  special  case  of  a  Gauss-Markov  process.  The  system  model  is  the  one 
of  Chapter  2  Section  2.  The  additional  Linear-Gaussian  assumption  yields, 
as  in  Chapter  5,  to  an  explicit  determination  of  the  distributions 
P(x(k),d  jz  )  and  P(x(k),d  [ z  )  of  the  last  chapter.  Furthermore,  under 
this  assumption,  pruning  rules  R1  and  R2  can  be  implemented  in  a  simple 
way.  Finally,  an  additional  deletion  technique,  specific  to  Gauss-Markov 
processes,  can  be  found. 


7.1  Computation  of  some  densities  under  a  Linear- 
Gaussian  assumption 

k  k 

Under  the  Linear-Gaussian  assumptions  (2.4)-(2.6),  P(x(k) ,d  jz  ) 

.  ki  k 

and  P(x(k),d  |z  ),  are  normal  functions,  weighted  by  some  coefficients. 
Here,  we  derive  how  these  functions  can  be  easily  computed  recursively. 
We  have  the  relation: 


P(x(.k)  ,dkj  zk)  »  PCxCk)jd\zk)  x  P(dk|zk) 


(7.1) 


k  k 

The  distribution  P(x(k)jd  ,z  )  represents  the  conditional  a  posteriori 
density  for  the  continuous  state  x(k)  on  a  discrete  state  sequence 

It 

d  «  (d(0),...,d(k)).  According  to  Section  3  of  Chapter  S,  we  know 

i  *k  k  k 

that  P(x(k)|d  ,z  )  is  a  Gaussian  density.  Its  mean  )  and  covari- 

k  k 

ance  P^j^Cd  )  depend  on  the  discrete  sequence  d  through  equations 

(2.4)-(2.6).  They  are  given  recursively  by  the  following  equations, 
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obtained  directly  from  the  Kalman  Filtering  Equations  (5.10) -(5 . 15) . 
Here,  the  state  transition,  observation,  and  covariance  matrices  at  time 
k  depend  also  on  the  discrete  state  values  d(k)  ■  p,  d(k+l)  *  q  at  time 
k  and  k+1: 

*k|k«k’  ’  %|k-lCd^  *  ^  (7.2) 


Vl|k(dk)  ■  fCP.I.^^IkCd11)  (7.3) 

Pk|kCdk;)  "  f1  *  ^cABCP.WJFkik.ifA  (7-4) 

k+1  k  T 

Pk*l|kCd  >  a  FCP.q»k)Pk|kCd^P  (p.q,k)  ♦  Q(p,q,k)  (7.5) 

KfcOft  *  pk|k.iCA  (P.k)[HCp,k)Pkjk^Cdk)HTCp,k)  + 

R(P.k)]"1  (7.6) 

xQ| _i Cd0^  a  0  and  P0|.1W0)  «  PQ  is  the  initial 

condition  (7 . 7) 


F(P><bk),  G(p,k),  Q(P,q.k),  R(p,k),  ?Q  are  defined  by  (2.4)-(2.6). 

The  coefficient  in  front  of  the  Gaussian  density  is  equal  to: 

P(dk|zk)  *  P(dk,zk)/P(zk)  where  PCz^)  is  only  a  normalization  term, 
k  k 

P(d  ,z  )  can  be  determined  recursively  from: 

P(d\zk)  -  P(z(k)|zk'1,dk)  x  P(dCk)!d(k-l))  x 

P(dk'1,zk"1).  (7.8) 

P(d(k) jd(k-l) )  is  known,  and,  as  we  saw  in  Chapter  5,  P(z(k) | zk_1,dk) 
is  the  probability  density  for  the  innovation  process,  conditioned  on 
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d‘ .  Therefore  we  have: 


P(z(k)| zk-1,dk)  -  N[z(k)  -  HCp,k)xk,k_1(dk),  B (dk) ]  (7.9) 


where 


B(dk)  «  H(p,k)Pk(k_1(dk)HT(p,k)  ♦  R(p,k) 


(7.10) 


k,  k 

Thus,  P(x(k),d  |z  )  can  be  easily  computed  recursively  using  (7.1)- 
(7.10). 

Similarly,  we  have: 


P(x(k),dk|zk)  *  P(x(k) [ dk,zk)  *  P(dk| zk) 


(7.11) 


l  ,k 


According  to  Section  5  of  Chapter  S,  P(x(k)|d  ,:  )  is  a  Gaussian 

k 

density,  conditioned  on  d  ,  weighted  by  a  coefficient,  and  whose  mean 

k  k 

and  covariance  are  x^^d  )  and  Pkjk(d  )  are  given  by  (7.2)-(7.7): 


P(x(k)|dk,zk) 


/  (2ff)nj  pkjk(dk^ 
y  (2ir)nk{ P  ^ k(dk) 


x  N(x(k)  - 


Mk 


CO, 


(7.12) 


k  k  k  k 

P  1  (d  )  is  t.’  e  covariance  of  the  whole  state  sequence  x  ,  conditioned 

k  k  I  k  k 

on  d  .  The  computation  of  P  1  (d  )  can  be  done  as  follows:  the  diagonal 
k  I  k  k  k 

::terms"  in  P  1  (d  )  are  equal  to  | ^ C d  ),  the  smoothed  covariances  de¬ 
fined  in  Chapter  S  by  (S .23) -(5 . 24) ;  the  off-diagonal  matrices  P.  .  . 

1+J  9  1 

corresponding  to  the  (i+j)**1  row  block  and  it^1  comumn  block  are  equal  to: 
Pi+jji  ■  F(d(j-1) ,d(j) , j-1)  x  ...  x  FCd(i)>dCi*l),i)Pi|k(dk) 


'(F(d(j-1)  ,d(j)  ,  j-1)  x  ...  x  F(d(i),dCi*l),i)f.,k 


'  ^i>j!k)55i|k 


(7.13) 
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Thus,  P(x(k) ,d^| can  be  obtained  from  (7. 2) -(7. 15) . 

These  functions  are  used  for  the  pruning  rules  of  Algorithm  1  and  2 
of  the  last  chapter.  The  LineaT -Gaus si an  assumption  yields  a  simple  im¬ 
plementation  of  some  of  these  rules.  This  is  the  subject  of  the  follow¬ 
ing  section. 

7.2  Actual  Implementation  of  R1  and  R2  under  the  Linear-Gauss ian 
Assumption 

Pruning  rules  R1  and  R2  of  the  last  chapter  are  based  on  functional 

k  k 

inequalities  between  the  above  distributions,  P(x(lt),d  |z  )  for  R1  and 
k  k 

P(*00,d  jz  )  for  R2.  It  is  straightforward  to  find  a  simple  test  equi¬ 
valent  to  the  functional  inequality  between  two  Gaussian  distributions , 
weighted  by  some  coefficients.  This  is  the  subject  of  the  following 
Lemma: 

Cx-x.^B^Cx-x.) 

Lemma:  Let  f^'x)  *  c^exp  -{ - = - i - }  and 

fXrX  yfB'^X-X.) 

f2(x)  *  c^exp  -C - - - - - =—}  be  two  functions  defined  over  a  state 

space  X,  such  that  0  <  and  0  <  B,.  Then,  for  all  x  G  X, 
fjCx)  <,  f2(x)  <->  (  0  <  Bj  <  B2 

(  (x1-x2)TCB2-B1)"1Cx1-x2)  <_  2  in 

Proof:  The  proof  is  left  to  Appendix  B. 

This  Lemma  yields  a  very  simple  implementation  of  pruning  rule  R1 
and  R2  in  the  Linear-Gaussian  case: 
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k  k 

Pruning  Rule  R1(L):  Let  and  d^  be  two  discrete  state  sequences  such 
thac  dl(k)  >  d2(k).  Let  ^(dj),  Pk|k(d‘j,  xk|k(d^),  PkjkCd^  be  the 
corresponding  means  and  covariances  as  defined  by  (7.2)-(7.7).  Let 

e<dJ-4’  ■  (\| k<d2>  -  *k|ktdl',T(Pk|kW2)  - 

k  -1  k  k 

Pk|ktd?5  ^|k^P  -  *k|k(dl»  *  2  in  ♦ 

'Pele.fdbt  2 


in 


rklkL°2J 

’kik(d5>i 


(7.14) 


If  <  Pic[k^d2^  311(1  ®*-dl,d2^  —  0>  i^en  dk  can  be  immediately  dis¬ 

carded  at  time  k,  without  increasing  the  probability  of  an  error  in  es¬ 


timating  d* 


Proof..  The  proof  comes  directly  from  pruning  rule  Rl,  the  expression  of 
k  k 

P(x(k),d  |z  )  in  the  LG  case  and  the  applications  of  the  Lemma. 

We  have  also  the  following  pruning  rule  R2(L),  obtained  directly 
from  R2: 


Pruning  Rule  R2(L): 
that  dj (k)  *  d2(k). 
be  the  corresponding 
Let  also: 


k  k 

Let  d^  and  d^  be  two  discrete  state  sequences  such 

*k|k<^’  Pk|k‘d!^  pMl'(dl)-  *k|k<d2>’pk|k'd2>’ 
means  and  covariances  as  defined  by  (7.2)-(7.13) . 


k  k 

Y(dJ,dp 


■  «k]k(dk)  -  *k!k<d?” 
W2>  -  *k|k(d^>  «  2  *" 


'<pk[k'd5>  -  pkikcdl”- 


P(dkUk) 

pcdSlsV 


Jin 


jPklk(dk) 

Pklk(dk)" 


C.15) 


r 
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I 


If  pk J ^ Cd and  YCd^.d^)  i  0,  then  d^  can  be  immediately  dis¬ 
carded  at  time  k,  without  increasing  the  probability  of  an  error  in  esti- 

K  K  K. 

mating  the  whole  hybrid  state  sequence  s  ■  (d  ,x  ). 

Proof:  Hie  proof  is  similar  to  the  one  of  R1(L)  and  is  omitted. 

]( 

(d2)  means  that  we  can  only  discard 

distributions  with  "small"  covariances.  Of  course,  a  necessary  condi - 

k  k  k 

tion  for  considering  the  deletion  of  d^  by  d^  is  that  d^  be  at  least 

ki  k  ki  k 

more  likely  at  time  k,  i.e.,  P(d^|2  )  <_  P Cd^ ( 2  )•  Thus,  if  we  want  to 

jc 

discard  a  given  sequence  dQ,  we  try  to  compare  it  only  with  other  more 

k 

likely  sequences  d  ,  whose  associated  Gaussian  distributions  have  a 
larger  covariance  (see  Fig.  7.1). 

Unfortunately,  no  such  simple  tests  can  be  found  for  the  generalized 
pruning  rules  GR1  and  GRZ,  which  are  certainly  more  powerful,  in  terms 
of  the  number  of  discrete  state  sequences,  that  could  be  deleted  at  a 
given  time.  For  these  rules,  the  conditions  would  have  to  be  verified 
numerically. 

So,  there  exists  a  tradeoff:  complexity  vs.  performance  in  our  al¬ 
gorithms.  The  need  for  an  easy  implementation  of  the  pruning  rules 
leads  to  the  use  of  Rl(l)  or  R2(L),  which  are  less  powerful  than  GR1  or 
GR2,  in  the  sense  that  more  state  sequences  will  remain  a',  time  K. 

Another  interesting  feature  of  the  LG  case  is  the  following:  the 
future  values  of  the  coefficients  of  comparison  8  ory,  between  2  state 
sequences  following  the  same  path  from  the  current  time,  can  be  deter¬ 
mined  at  this  time,  independently  of  the  future  measurements.  This 
leads  to  some  other  sequence  deletion  opportunities ,  as  we  shall  next 


The  condition  P 


k  |  k 


(.<§ 


<  P 


k|k 


see. 
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FI  OURE  7.1 


pfxfiH'b*) 
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.3  Another  Pruning  Rule  Specific  to  the  LG  Case 


k  k 

Let  us  consider  two  state  sequences  and  d^  such  that  d^(k)  = 


k  K 

d2(k)  and  P^j^Cdp  <  P^j^Cd^.  If,  according  to  our  objective. 


Ir  Ir  k  k  k 

BCd^.dj)  >  9  or  YCd^.d^)  >  0,  then  we  cannot  discard  immediately  d^  at 


time  k,  using  pruning  rule  R1  or  R2.  However,  the  future  values 


8(d^  ,d2  )  or  Y(d^  ,d2  )  of  the  coefficients  at  time  k'  >  k,  where 


d*’  *  (d*,d(k+l),...,d(k'))  and  d!j\«  (d*  ,d(k+l) , . . .  ,d(k')) ,  are  pre- 


coraputable  at  time  k.  This  is  true  because  of  the  following  theorem. 


k  k 

Theorem:  Let  d^  and  d2  be  two  discrete  state  sequences  such  that 
d^(k)  a  d7(k)  and  let  B(d^,d2)  and  YCd^.dj)  be  the  corresponding  coef¬ 


ficients  as  defined  by  (7.14)-(7.1S) .  Assuming  F(p,q,k)  is  invertible, 
then  for  all  2(k+l) ,  for  all  d(k+l) , 

|B(d*+1)| 


•  ecdj.d^  ♦  An 


*1'.  tn 


B(d~  *) 


W*iCdi  3 


An 


"klk^ 


(7.16) 


k+i  k+i  k  k  |B(d^+1) 

Y(d*+1,d*+1)  -  Y(d*,d*)  *  An  - ♦ 


|pk+l|k+l(dk+l5 
An  — .  rrr— ; — rr? - An 


lnwf1)] 


[pkik^kj 


P 


k+li  k+l^k+1^ 


(7.17) 


Proof:  The  proof  of  this  theorem  is  in  Appendix  C. 


The  main  consequence  of  this  theorem  is  that  we  can  predict,  through 


the  coefficients  0  or  y.  the  future  evolution  in  the  relative  position  of 
two  distributions ,  associated  with  a  common  sequence  from  the  current 
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time,  since  the  covariances  are  pre computable  at  time  k. 

Furthermore,  using  (7.15)  and  (7.16),  we  have  the  following  inequali 

ties : 


eCdJ+l,d!|*1)  <  0(d*,d£)  (7.18) 

Y(d*+1,d*+1)  <  Y(d*,d*)  (7.19) 

lc  lc  V  ^1  V+l  IT 

So,  starting  at  time  k,  the  sequences  Btd^.dj).  B(d^  .d*  ) , . . .S(d^ ,d^) 

or  Y(dk,dk)  >  YCd^.d^1 ),..., YCdJ.d^),  where  for  all  i  G  [k,K], 

dd  »  (dk,d(k+l) , . . . ,d(i))  and  d*  *  (dk,d(k+l)  , . . .  ,d(i)),  are  decreasing 

jj 

sequences.  Even  if  we  cannot  discard  d^  at  time  k  using  R1(L)  or  R2(L) 

because  8(dj,d2)  >  0  or  Y^.dj)  >  0,  there  is  a  chance  that  B(dk  ,d^  ) 
k '  k' 

or  Y(d^  ,d2  )  become  negative  at  a  later  time  k'  that  we  are  able  to 

k' 

determine  through  (7.16)-(7.17) .  So,  we  can  delete  d^  immediately  at 
time  k,  if  the  two  following  conditions  are  satisfied: 


Ci)  PkikCdi5  <.Pk|k(d2^  (which  implies  Pk,|k, (<£')< 


or 


.  ,  k'  |B(dJ')J 

(ii)  0  (dk  dk)  <.  Z  In- - i- 

1  2  j-k+1  l  B  (dp  | 

X 


B(dh| 


iPk-!k'(di'}1 
♦  In  *  l  --p —  ♦ 


' k*  k 


,cd;  )i 


k »  k  ’  k 1 
PK  |K  (d?  ) 


v  If  *  I  J  |  I  •  (.“n  ;  I 

Y(d  d,)  <_  Z  - — rr  *  &n  i  ftf — FT* 

1  2  j-k*l  i B(d^) i  !p  1  (d,  )i 


Jin 


[pkik(d2) 

^7 


(7.20) 


(7.21) 
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Furthermore,  we  can  simply  eliminate  the  state  sequence  d^  at  time 
k,  if  for  all  the  future  sequences  d^(k+l) , . . .  .d^CK) ,  we  can  find  at  least 
one  sequence  dk  *  (d^d^k+l) , . . .  .d^K))  such  that  PklkCdk)<  Pk)k(dk) 
and  (7.20)  or  (7.21)  is  satisfied  for  k'  *  K. 

This  constitutes  the  additional  pruning  rule  in  our  graph.  Let  us 
state  it  formally  for  both  objectives  as  follows: 


If 

Additional  Pruning  Rule  AR1(L):  Let  d^  be  a  given  state  sequence  up  to 

£ 

time  k.  If,  for  all  dg(k*l) , . • . ,dg(K) ,  there  exists  a  sequence  d  = 
(d\d0(k*l),...,d0(K))  where  d(k)  -  dQ(k)  and  Pk[k(dk)  <  P  j  (dk)  and 

[W# 


Tfa.Jh<  j  *JB(dn] 


j=k+l 


|B(dJ)| 


+  An 


K[  K 


.tn  |pk|k'dk) 


(7.22) 


Then  dQ  can  be  immediately  deleted  at  time  k,  without  increasing  the 
probability  of  an  error  in  estimating  dK. 

We  have  the  following  pruning  rule  for  our  second  objective: 


Additional  Pruning  Rule  AR2(L):  Let  d^  be  a  given  state  sequence  up  to 

K 

time  k.  If,  for  all  dg(k+l) , . . . ,dg(K) ,  there  exists  a  sequence  d  » 
(dk,d0(k+l) . d0(K))  where  d(k)  -  dQ(k)  and  Pk|k(dk)  <  Pk|kCdk)  and 

k  k  £  1*C4)1  |PK,K(dj)|  |p*l*fA| 

0  “j-fel  i B(d^)  j  !PK!K(dK)|  |P  1  (dj)l 

(7.23) 


V 

then  d^  can  be  immediately  discarded  at  time  k,  without  increasing  the 
probability  of  an  error  in  estimating  the  hybrid  sequence  s  . 

When  the  covariances  tend  quite  rapidly  (in  a  few  steps)  to  their 
steady-state  value,  the  coefficients  S  and  y  tend  also  rapidly  to  a 
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steady-state  value  on  a  given  path.  This  means,  that  in  this  case  these 
coefficients  must  be  positive  but  close  to  0,  in  order  *ot  AR1 (L)  or 
AR2(L)  to  work  effectively. 

We  illustrate  AR1(L)  in  figure  (7.2),  for  the  case:  d(k)  can  take 
only  2  values  { 0 , l)  at  each  time  k  and  k  *  K-l. 

Together  with  R1(L)  and  R2(L),  these  pruning  rules  lead  to  the  fol¬ 
lowing  algorithms  that  are  actually  implementable  versions  of  Algorithm 
1  and  2  of  the  last  chapter. 

7.4  State  Estimation  Algorithms  under  the  Linear-Gauss ian  Assumption 

To  apply  R1(L)  or  R2(L) ,  we  need  to  compute,  at  each  time  k,  the 
probability  of  each  remaining  discrete  sequence  at  that  time,  its  asso¬ 
ciated  continuous  state  estimate  and  covariance;  for  AR1(L)  and  AR2CL), 
we  need  also  to  compute  the  future  covariances . 

At  time  K,  we  must  choose  the  discrete  state  sequence  d  among  all 
the  remaining  sequences.  Then,  for  the  hybrid  state  estimation  algo- 

K 

rithm  2(L),  an  estimation  of  the  continuous  state  sequence  x  can  be 
simply  done,  using  the  Kalman  Smoothing  Equations  of  Giapter  5,  con¬ 
ditioned  on  d^  *  (d(0) , . . .  ,d(K))  ;  the  estimate  x*^  ^(d1^  *  (x^j  ^Cd*S  , . . . , 

$K|K(d  ))  is  given  recursively  by  the  following  equations,  where 
d|(k)  *  p,  d(k+l)  »  q:  for  k  ■  K-l.K-2, . . .  ,0, 

XKJ K(d  )  is  the  boundary  condition  for  k  ■  K-l 

*  ^k^cp.q.k^UkCa11*1) 


(7. 25) 
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WJK>  •  pk|ktak>  *  <?•“> 

*k|k'^'  *k.l|kWk“)'  Pk|k(ak)'  ''k.llkt?*1)  *r5  S1™'”'  C7-2)-(7.TJ. 
Using  all  the  above  results,  we  have  the  following  alforithm  for 

K 

estimating  the  discrete  state  sequence  d  ,  according  to  the  first  oojec- 
tive  of  Chapter  6: 


ALGORITHM  1(1) 


STEP  0: 


STEP  k: 


*k-i|k-i«k'l>-  pk.I|k.i'dk'15 

P(dk“1|zk'1)  ,  for  all  dk-1  GS^1  . 


RECEIVE  z£k) 

COMPUTE  Xj^j  k(dk) ,  Pk|k(dk),  P(dk|zk)  for  all  (£  G  Q^”1  x  d. 
APPLY  R1(L):  for  all  dk  G^'1  x  d,  there  exists 
dk  Q&P'  ^  x  d  .  {dj^}  such  that: 

ki  k.  ki  k 

a)  pc«ln  i  pcdK{ 

NO 


YES:  Continue 
NO:  Store  dk  in  3k 


(ii) 


pkik'd5’ 


YES 

< 

NO 


pkik'dk> 


YES:  Continue 
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NO:  Store  dk  in 
YES 

(iii)S(dk  dk)  <_  o 
u  NO 

YES:  Delete  dk 

NO:  Continue  to  AR1(L) 

APPLY  AR1(L):  for  all  dQ  Ck-^l) ,  •  •  •  ,dQ(K)  ,  there  exists  dk  verifying 
(7.22) 

YES:  Delete  dk 
NO:  Store  d£  in  3  s 

k  ■  k+1  and  repeat  until  k  *  K. 

STEP  K:  (2IK 

(  P(dK|  z*S  for  all  dK  S  0K  . 

Find  Sup  PCd^t1^,  The  corresponding  state  sequence  is 
■*K 

equal  to  d  .  END. 

We  have  also  Algorithm  2(L)  for  estimating  the  hybrid  state  sequence 
sK  *  (dK,xK) : 


ALGORITHM  2  CL) 


■  do 

Vi  * 
y  poi-i  ■  po 


STEP  0: 
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STEP  k: 


1 

)  $  fJ<-l<,  p  fJt-l..  0k-l  k-1  k-1. 

\  \-l\k-l<-d  )p  Pk-l|k-lCd  hP  (d  )p 

'  PCd*'1^'1)  for  all  dk-1  6  S2k_1 


RECEIVE  z(k) 

COMPUTE:  %jkCdk),  Pk[k(dk),  Pk'k(dk),  P(dkj sk)  for  all 
dk  G^"1  x  D. 


APPLY  R2(L):  for  all  dk  6  0k_1  *  D,  there  exists  dk  G^'1  x 
D  -  {dk}  such  that: 

(i)  P(dJ|2k/lS  P(dk|2k) 

U  NO 


YES:  Continue 
NO:  Store  dk  in  2^ 
YES 

Cl»  pk|k«5>  i  pk|k(dk) 


YES:  Continue 

NO :  Store  d^  in  S* 

.  ,  YES 

(iii)  Y(d*dK)  <  0 

u  NO 

YES:  Delete  djj 

NO:  Continue  to  AR2(L) 


APPLY  AR2(L):  for  all  dQ(k+l) * . • .,d0Qc) ,  there  exists  dk  veri¬ 
fying  (7.23) 

YES:  Delete  dk 
NO:  Store  dk  in  0k 
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k  =  k+1  and  repeat  until  k  =  K. 

STEP  K: 

|  PKI  V),  P(dK|zK)  for  an  dK  E3K 
(  XK|  KCd*S 

Find  Sup  , ,  i.  „  x  PCd^lz^).  The  corresponding  discrete 

d*-S2r  VI  ?  COl 

state  sequence  is  equal  to  d^.  Finally,  «  (d^,  X^'^d^) 
where  K(3*S  *  {^j^Cd^),  k  =  K-l,  K-2,...,0)  is  given  re¬ 
cursively  by  (7. 24) -(7. 25) .  END. 

Algorithm  1(L)  and  2(L)  are  easily  implementable.  We  see  that 

£ 

both  objectives  do  not  necessarily  yield  the  same  estimate  of  d  , 

K 

since,  in  the  first  case,  we  choose  the  most  likely  sequence  d  ,  and 

v 

in  the  second- one,  we  pick  the  sequence  d  for  which 
maximum,  which  is  not  necessarily  the  most  likely. 


and  Conclusion 


Under  the  Linear -Gaussian  assumption,  pruning  rules  R1 (L)  and 
R2(L)  used  for  estimating  the  discrete  or  hybrid  state  sequence  of  a 
Markov  process  are  easily  implementable.  Another  pruning  rule  specific 
to  the  LG  case  has  also  been  derived.  These  implementable  pruning  rules 
lead  to  the  derivation  of  algorithms  1(L)  and  2(L),  which  are  more  re¬ 
stricted  than  Algorithm  1  and  2,  but  are  actually  implementable. 

A  further  problem  would  be  to  compare  the  number  of  sequences  we 
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could  delete  optimally  this  way,  with  sufaoptimal  algorithms  that  simply 
drop  very  unlikely  sequences  during  the  graph  formation. 

In  the  next  chapter,  we  apply  algorithm  1(L)  to  a  simple  failure 


detection  problem. 
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CHAPTER  3 

EXAMPLE:  APPLICATION  TO  FAILURE  DETECTION 


The  purpose  of  this  chapter  is  to  apply  Algorithm  A1C.L)  to  a  simple 
failure  detection  problem.  In  particular,  we  will  look  at  the  efficiency 
of  our  optimal  pruning  rule  R1C.L). 

8.1  A  Failure  Detection  Problem 


We  consider  a  one-dimensional  dynamic  system,  whose  continuous  state 
is  observed  through  a  sensor  subject  to  intermittent  failures  occurring 
at  random  times:  in  the  normal  case,  the  sensor  provides  state  observa¬ 
tions;  in  case  of  failure,  only  a  burst  of  noise  is  received.  The  role 
of  the  failure  detection  system  is  to  decide  upon  the  time  and  duration 
of  sensor  failure,  with  a  minimum  probability  of  an  error:  this  is  an 
alarm  task. 

We  apply  the  single  hybrid  Gauss-Markov  model  of  Chapter  2  Section 
2  to  this  problem,  and  pick  some  numerical  values. 

The  discrete  state  d(30  of  the  system  at  time  k  can  take  two  values: 

d(k)  ■  ( 0  if  the  sensor  has  not  failed 
(l  if  the  sensor  has  failed. 

The  continuous  state  x(k)  belongs  to  R.  We  have: 

-  the  discrete  dynamics: 

P(d(k+1)  *  lid(k)  ■  0)  ■  Q.l 

P(d(k+1)  *  0| d(k}  «  1)  ■  0.1 


0.1 


-  the  continuous  dynamics:  for  d(k)  *  0,l,x(k+l)  »  x(k)  *  w(k) 
where  w£k]  -  N(0,1) 

-  the  measurement  equation: 

if  d(k)  a  0,  then  2(k)  *  *00  ♦  VQ00  where  vQCk)  ~  N(0,1) 
if  d(lc)  *  1,  then  zQO  =  v^}  where  v^Ck)  -  N(0,50) 

Since  the  task  of  the  failure  detection  system  is  only  to  detect 
the  sensor  failures,  we  use  algorithm  A1(L1,  which  estimates  the  discrete 
state  sequence  of  our  system. 

8.2  Results 

The  system  is  observed  during  an  interval  [0,8].  We  assume  that 
a  noise  burst,  resulting  from  a  sensor  failure,  arises  at  k  =  2,5,4. 

We  study,  for  this  case,  the  efficiency  of  priming  rule  R1(L]. 

We  take  the  following  values  for  the  sequence  of  observations: 


k 

0  12  3  4 

5 

6 

7  8 

zOO 

3  2.5  15  -10  10 

0 

3 

2.5  2.5 

The  application  of  pruning  rule  R1(I0  in  the  gTaph  of  the  possible 
discrete  state  sequences  yields  the  diagram  shown  in  Fig.  8.1 

We  see  that  the  algorithm  is  working  pretty  well.  At  time  3,  we 
are  left  with  only  11  different  discrete  state  sequences,  using  R1(L) 


ilurb 
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(the  maximum  number  of  possible  sequences  is  512).  Furthermore,  the 
algorithm  confirmed  at  time  4  and  5,  that  a  failure  has  occurred  at  time 
2  and  3  respecitvely.  Therefore,  it  responds  quickly  to  the  sensor 
failure. 

So,  the  application  of  A101)  to  the  multiple  hypothesis  filter- 
detectors,  mentioned  in  Chapter  3,  should  lead  to  an  optimal  procedure 
for  detecting  a  failure,  that  keeps  the  amount  of  computation  at  a 
reasonable  level. 

This  example  highlights  the  efficienty  of  our  pruning  rule  in  a  very 
simple  case.  Of  course,  more  complete  simulations  would  have  to  be  run, 
to  estimate  the  real  capability  of  our  algorithm. 
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CHAPTER  9 

STATE  ESTIMATION  OF  A  MULTIPLE  HYBRID  MARKOV  PROCESS 


In  part  II,  the  state  estimation  of  a  single  hybrid  Markov  process 
has  been  studied.  We  now  consider  the  more  complex  problem  of  estima¬ 
ting  the  state  sequence  of  multiple  hybrid  Markov  processes  in  parallel. 

9.1  Problem  Statement 

The  multiple  system  model  considered  here  is  the  one  described  in 
Chapter  2  Section  3. 

As  for  the  single  system  problem  of  Chapter  6,  there  are  two  pos¬ 
sible  objectives  for  a  state  estimation  of  a  multiple  hybrid  system, 

according  to  the  application.  One  possible  objective  is  the  MAP  esti- 

K  K  1C 

nation  of  the  total  hybrid  state  sequence  s  •  C<i  ]  of  the  multiple 
system.  In  other  cases,  e.g.,  in  some  multi-target  tracking,  we  could 
be  more  interested  in  the  number  of  targets  in  a  certain  area,  and  thus 

K 

would  use  an  MAP  estimation  of  the  discrete  state  sequence  d  of  the 

K 

multiple  systems.  Then,  conditioned  on  the  estimate  of  d  ,  we  can  look 

(; 

for  the  MAP  estimate  of  the  continuous  state  sequence  x  . 

Because  of  the  Markov  and  memoryless  properties  C2.9)-(2.11) ,  we 
can  apply  directly  algorithm  1  or  2  of  Chapter  6  to  the  multiple  system, 
which  is  the  cartesian  product  of  the  individual  systems.  However,  such 
an  algorithm  would  not  take  advantage  of  a  possible  decomposition  of  the 
problem,  due  to  the  relative  independence  among  the  subsystems.  If  we 
could  find,  in  both  cases,  the  optimal  estimates  by  running  in  parallel 


p  independent  algorithms  1  or  2,  one  for  each  subsystem,  then  the  amount 
of  computation  would  be  gTeatly  reduced:  the  total  number  of  possible 
sequences  in  the  corresponding  p  independent  graphs  is  p  *  N  ,  which  is 
much  smaller  than  the  number  of  possible  sequences  in  the  corresponding 
graph  for  the  multiple  system,  equal  to  (N1*)*. 


9.2  Decomposition  of  the  Estimation  Problem 

The  purpose  of  this  section  is  to  find  the  assumptions  that  would 
make  the  general  objective  of  the  estimation  problem  equivalent  to  a 
set  of  local  independent  objectives,  in  order  to  get  a  decomposition 
of  our  algorithm. 

As  mentioned  in  Chapter  2  Section  3,  the  MAP  hybrid  state  estimate 
of  sK  given  by  sK  for  which  P(sk|  Y*) ,  or  equivalently  P(YK|sK)  *  P(sK) 
is  maximum.  Using  the  Markov  and  memoryless  properties,  we  have: 

K 

PCYK|s*S  X  PCs*)  »  JI  P(YCk)|sCk))  x  P(sCk)lsCk-l))  (9.1) 
k*0 

Similarly,  for  the  discrete  state  estimation  of  the  multiple  system, 

AK  A  A 

the  objective  is  to  find  d  *  (d(0) , .. .d(K))  for  which: 

P(YK!dK)  x  P(d *)  .  JI  PCYCk)idk,Yk_1)  x  p(d(k)  |  d(k-l))  (9.2) 
k-0 

is  maximum. 

We  recall  the  independent  dynamics  and  independent  measurement 
generation  processes  assumptions  that  we  made  in  Chapter  2: 

Assumption  (A.l):  The  dynamics  of  the  subsystems  are  assumed  to  be 
independent,  i.e.,  for  all  k  ■*  0,1,..., K, 


P(sCk)jsCk-l) 


(9.3) 


P 

n 

i=0 


PCs.OOis^k-1)) 


Assumption  (A. 2):  The  sets  of  measurements  generated  by  the  subsystems 
and  the  sets  of  false  alarms  are  assumed  to  be  statistically  independent, 
i.e.,  for  all  k  »  0,1,. ...K  , 


P(Y1(k),...,Yp(k),Yf(k)!s(k)) 


P 

n  P(Y  (k)k(k))  x  P(Y  (k)) 
i*l 


(9.4) 


This  implies,  that  for  all  k  ■  0,1,..., K, 


P(YX (k) , . . . ,Yp(k) ,Y£(k)  [  dk,Y*_1, . . . ,Yp_1 .Y**1)  * 

I  P(Y (k)|dk,Yk_1)  x  P{T  (k)  )  (9.5) 

i-1  1  1  1  r 

At  this  point,  if  the  partition  of  Y(k)  into  { Y  (k) , . . . ,Y  (k) , 

i  p 

Yf(k)}  were  a  priori  known,  then,  under  assumptions  (A.1)-(A.2),  our 
estimation  problems  would  be  equivalent  to  find  for  each  i  ■  l,2,...,p, 

r  K 

the  sequence  s:  for  which  II  P(Y.  (k)|  s.  (k))  x  P(s.  (k)js.(k-l))  is 

1  k»0  1  1  1  x 

maximum,  or  under  the  second  objective,  the  sequence  d^  for  which 
K  k  k-i 

H  P(Y, (k) j d. ,Y?  )  x  P(d. (k) | d. (k-1))  is  maximum.  The  estimation 

k*0  x  11  1  1 

problem  would  be  solved  by  p  independent  subproblems. 

Unfortunately,  as  we  emphasized  in  Chapter  2,  in  general  this  par¬ 
tition  is  unknown.  We  can  only  make  a  hypothesis  for  the  partition  of 
Y(k)  at  time  k,  that  is  called  a  data  hypothesis  at  time  k,  denoted  by 
A(k) .  We  saw  that  the  probability  distribution  for  the  set  of  measure¬ 
ments  Y(k)  is  equal  to: 
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P 

P (Y (k) !  s 00 )  s  I  H  P(T.(A£K))|  s.(k))  X  PCY  (A(k))  * 

{ A(k)}  i-1  1  1  f 

P£A(k)jd(k))  (9.6) 

and 

PCY(k)|dk,Yk'1)  »  Z  I  P(Y.(A(k))|dk,  Yj^CACk)))  x 
{  Ak>  i-1  1  11 

P(Yf(A(k))  x  P(Ak|dk)  ,  (9.7) 

where  Y^Afk)) , . . . ,Yp(A(k)) ,  Y^(A(k))  denotes  the  hypothesized  parti¬ 
tion  of  Y(k)  under  A(k). 

Thus,  even  though  the  system  dyianics  and  measurement  processes 
are  independent,  the  data  association  uncertainty  introduces  a  depen¬ 
dency  among  the  subsystems  for  the  estimation  problem.  If  we  cannot  re¬ 
duce  the  uncertainty  in  the  data  hypotheses  by  some  further  assumptions, 
then  no  decomposition  is  possible. 

As  mentioned  above,  an  a  priori  knowledge  of  the  partition 
Yj (k) , . . . ,Yp(k) ,Y^(k)  would  lead  to  this  decomposition,  since  it  would 
imply  that  only  one  data  association  is  possible  and  (9.6)-(9.7)  would 
reduce  to  (9. 4) -(9. 5)  respectively. 

Some  weaker  assumptions  can  also  lead  to  the  decomposition.  For 
this  purpose,  we  define  the  following  regions  of  the  measurement  space 
Z. 

Definition:  A  measurement  region  R^(k)  is  defined  as  the  set  of  mea¬ 
surements  that  can  a  priori  come  from  subsystem  i  at  time  k. 

In  some  digital  estimation  problems,  it  can  correspond  to  the  pos¬ 
sible  discrete  values  of  the  measurement  emitted  by  a  subsystem.  In 
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multitarget  tracking  problems,  as  we  shall  see  in  Chapter  11,  it  cor¬ 
responds  to  a  region,  outside  which  no  measurement  can  arise  from  the 
target . 

If  at  each  time  k,  the  sets  YCk)  A  R^(k)  are  mutually  exclusive, 
then  the  data  association  uncertainty  is  greatly  reduced.  A  partial 
partition  of  the  set  of  observations  Y(k)  is  possible:  the  measurements 
located  outside  any  regions  R^OO  can  definitely  be  assigned  to  the  set 
of  false  alarms.  Each  measurement  inside  R^(k),  however,  can  be  assigned 
either  to  subsystem  i  or  to  the  set  of  false  alarms.  If  the  detection 
and  false  alarm  processes  are  such  that  the  numbers  of  detections,  missed 
detections,  and  false  alarms  in  a  given  region  are  statistically  inde¬ 
pendent  of  the  numbers  in  any  disjoint  region,  then  nothing  correlates 
anymore  the  different  subsystems  and  the  decomposition  of  our  estimation 
problem  should  follow.  This  assumption  is  weaker  than  the  one  that  would 
suppose  an  a  priori  knowledge  of  the  sets  Y^  00  , • • • .Y^OO  ,  Yf (k) ,  since 
there  remains  a  data  association  uncertainty  inside  each  region  R^(k), 
between  actual  observations  and  false  alarms.  CSee  Fig.  9.1) 

This  independence  assumption  can  be  stated  formally  as  follows. 

If  at  time  k,  the  sets  Y(k)  H  R^OO  are  mutually  exclusive,  then  we 
partition  Y(k)  in: 

P 

YGO  -  {Y(k)AR  00.....Y00  H  R  00.  Y(k)  f»  [Z  -  u  RiOOJ}. 

1  P  i»l 

P 

A  measurement  in  Yfk)  fi  [Z  -  u  R. (X)]  is  associated  to  the  set  of 

i-1  1 

false  alarms  Y  (k) ,  since  it  lies  outside  any  region  R.  00  •  In  this 
f  i 

case,  the  data  hypothesis  A(k)  takes  the  form:  A(k)  *  (A^Ck)  , . . .  ,Ap(k) , 
A^OO),  where  Ai(k)  is  thrt  corresponding  data  hypothesis  from  d^k) 
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into  Y(k)  A  R. (k) ,  and  A. 00,  and  A.(k)  is  the  association  of  the  mea- 

1  t  j  r 

surements  in  Y(k}  A  (Z  -  U  R.00)  to  the  set  of  false  alarms  Y  00  • 

i*l  1  £ 

We  make  the  following  assumption: 

Assumption  (A.  31:  If,  at  time  k,  the  sets  Y(k)  A  R^0O  are  mutually 

exclusive,  then  we  assume  that  all  the  data  hypotheses  A.  Q0,...,A  00, 

P  P 

A_00  are  statistically  independent,  i.e.,  P(A(k)|d(k))  =  II  P(A. (k) | 

i*l  1 

d.(k))  x  P(Af0O). 

In  a  multitarget  tracking  application,  if  the  detectability  of 
each  target  is  independent  of  the  others,  and  if  the  process  generating 
a  false  alarm  is  independent  of  the  other  false  alarms,  then  assump¬ 
tion  (A. 3)  is  satisfied.  But,  as  it  is  often  the  case,  if  the  number 
of  false  alarms  in  the  area  is  assumed  to  follow  a  Poisson  distribution 
then  Assumption  (A. 3]  would  not  ...Id;  such  a  model  for  the  false  alarms 
would  prevent  any  decomposition  of  our  problem. 

Thus,  if  the  regions  Y(k)  A  R^(k)  are  mutually  exclusive,  then 
the  possible  data  hypotheses  take  the  fora:  A(k)  *  {A^0t),...,A  (k), 
A^fk)}  defined  above.  In  this  case,  we  have  for  the  hybrid  state  es¬ 
timation  problem  under  (A. 2) -(A. 3), 

P(Y(k)|s(k))  -  Z  P CY Ck) 1 s Ck) , A (k) )  x  P(A(k) | dCk)) 

(A(k)} 

P 

«  s  n  pcyoo  n  R.(k)|s.ck),A.oo)  x 

{A1(k)>,...,{ApW}i-l  ill 

P(A.(k)|d.{k))  x  P(ZfCk))  x  PCAf(k))  *  PCZfCk))  x 
P 

P(Af(k))x  B  Z  P(Y(k)  A  R.(k)|s.(k),A.Ck)) 
f  i-l{AiCk)}  1  1  1 


X 
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P(A.Ck)|d.(k))  =  P(Z£00  x  P(AfCk))*II 

i=l 


PCY(k)n  RiCk)|si(k)) 


p 

where  Zf(k)  denotes  Y(k)  n  [Z  -  U  R.(k)]. 

*  1 

So,  using  (A.l),  the  problem  of  finding  sK  for  which  P(sK|y^  is 
maximum,  is  strictly  equivalent  to  the  p  following  independent  sub- 
problems:  for  all  i  »  l,2,...,p,  find  s£  for  which: 


K 

H  PCY(k)  r\  R.  (k)|s.(k))  x  PCs.  (k)is.Ck-l))  is  maximum;  and 
k=0  11  11 


Similarly,  for  the  discrete  state  estimation  problem,  we  can  show 
under  the  same  assumptions  that: 


P(Y(k)|  Y^.d*)  - 


PCZfCk))  *  PfAf Ck) )  x 


p 

n  p(Y(k)  n  r.  oo | 

i»l  1 


Therefore,  the  problem  of  finding  d^  for  which  P(d^|Y*S  is  maximum  is 

equivalent  to  the  following  independent  subproblems:  for  all  i  *  1,2,..., 
•'K 

p,  find  d^  for  which: 

K  ,  . 

n  P(Y(k)  n  R. Ck) j  Y.  ,dv)  x  P(d. (k) j d. (k-1))  is  maximum; 
k*0  ill  li 

and  dK  * 

We  summarize  the  above  results  as  follows: 


Main  Result:  Under  assumptions  (A. 1)- (A. 2) -(A. 3) ,  if,  at  each  time 
k  *  0,1,..., K,  the  sets  Y(k)  O  R^{k)  are  mutually  exclusive,  then  the 
multiple  hybrid  and  discrete  estimation  problems  are  equivalent  to  a 
set  of  p  independent  subproblems,  one  for  each  subsystem. 


4 

Ai 
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Thus,  instead  of  running  one  algorithm  for  the  multiple  system, 
we  would  run  p  algorithms  in  parallel,  similar  to  the  ones  of  Chapter 
6.  This  would  lead  to  an  enormous  simplification,  since,  before  any 
deletions,  the  total  number  of  possible  state  sequences  at  time  K  would 
be  reduced  from  to  p  *  N^. 

If  some  of  the  sets  Y(k)  A  R^Ck)  do  overlap  during  the  period  [0,K], 
then  a  partial  decomposition  of  the  problem  is  still  possible,  by  run¬ 
ning,  as  in  Reid's  algorithm  [9],  one  algorithm  for  the  cluster  of  the 
corresponding  subsystems ,  and  independent  algorithms  for  the  remaining 
subsystems . 

A  further  problem  is  to  know,  if,  after  the  overlapping  of  some 
sets,  we  can  still  recover  the  decomposition  property  of  our  algorithm, 
when  they  separate  again.  This  is  the  subject  of  the  following  section. 

9.3  Actual  Implementation  of  an  Algorithm 

An  algorithm  for  estimating  the  hybrid  or  discrete  state  sequence 
of  a  multiple  Markov  process  must  take  advantage  of  any  decomposition 
opportunity  for  reducing  the  amount  of  computation. 

The  multiple  system  is  observed  during  an  interval  [0,KJ.  As  long 
as  no  overlap  of  any  sets  YCk)  A  R^k)  has  ever  occurred,  we  can  run  p 
algorithms  in  parallel.  If,  at  a  given  time  k,  some  overlapping  takes 
place,  then  we  can  still  revert  to  independent  algorithms,  if,  using 
the  pruning  rules,  we  find  a  single  remaining  survivor  through  the 
interval  of  overlap,  since  this  confirms  the  multiple  system  state 
during  this- interval .  Otherwise,  the  resulting  data  association  un¬ 
certainty  introduces  a  correlation  among  our  subsystems,  that  is  going 
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to  last  until  time  K.  So,  after  that  time  k,  one  algorithm  has  to  be 
run  for  the  overall  system. 

In  many  cases,  the  number  of  remaining  state  sequences  at  time  K 
is  still  going  to  be  well  beyond  our  storage  capabilities.  Then,  the 
algorithms  are  truncated  to  some  manageable  length  K^,  much  smaller 
than  K,  in  order  to  keep  the  number  of  hypothesized  sequences  at  a 
reasonable  level.  In  these  K^-scan  algorithms,  for  each  time  k  »  K^, 
Kr+1,...,K,  we  compute  recursively  the  discrete  or  hybrid  state  estimate 
at  time  k-K^  as  follows:  having  decided  upon  the  state  estimates 
d[k-Kr-l]  or  Sfk-K^-1]  at  time  k-1,  we  determine  the  most  likely  dis¬ 
crete  or  hybrid  state  sequence  from  k.- 1  to  k,  using  the  pruning 
rules  of  Algorithm  1  or  2  of  Chapter  6.  Then,  the  corresponding  dis¬ 
crete  or  hybrid  state  dfk-K^]  or  s£k-Kr]  on  this  sequence  is  chosen  as 

the  state  estimate  at  time  k-K  .  Since  at  time  k-1,  the  state  at  time 

r 

k-Kr-l  has  been  confirmed,  we  can  run  p  independent  algorithms,  if, 

from  to  k-K  to  k,  the  observation  sets  defined  in  Section  2  are 
r 

mutually  exclusive.  Thus,  this  guarantees  a  decomposition,  even  though 
there  has  been  some  overlap  before  time  k-Kp.  (See  Fig.  9.2). 

This  kind  of  situation  often  arises  In  multitarget  tracking.  When 
the  surveillance  system  has  just  been  initialized,  we  have  a  very  poor  pic¬ 
ture  of  the  area  with  much  uncertainty  in  the  location  of  the  dif¬ 
ferent  targets.  After  a  while,  the  target  state  estimates  become  more 
accurate  and,  if  they  are  largely  spaced,  the  regions  Y(k)  D  I^Ck)  be¬ 
come  mutually  exclusive.  Then,  some  crossings  can  occur  which  intro¬ 
duce  again  a  correlation  between  the  different  target  state  estimations 
for  a  non-zero,  but  finite,  duration. 


J 
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9.4  Summary 

The  computational  complexity  of  an  algorithm  for  estimating  the 
discrete  or  hybrid  state  sequence  of  a  multiple  Markov  process  led  us 
to  consider  a  decomposition  of  this  problem  into  several  independent 
subproblems.  In  order  to  get  this  decomposition,  some  model  assump¬ 
tions,  as  independent  dynamics,  independent  measurements,  independent 
detection  and  false  alarms,  had  to  be  made,  and  a  condition  on  the 
sets  Y(k)  n  R^k)  had  to  be  met.  For  the  application  to  multitarget 
tracking  in  Chapter  11,  the  actual  computation  of  these  regions  R^Ck) 
will  have  to  be  done.  This  will  show  that  our  multitarget  tracking 
problem  under  the  above  assumptions  can  be  equivalent  to  p  independent 
single  target  tracking  problems. 

We  next  consider  an  additional  linear-Gaussian  structure  on  the 


multiple  Markov  process. 
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CHAPTER  10 

STATE  ESTIMATION  OF  A  MULTIPLE  GAUSS-MARKOV  PROCESS 

This  chapter  is  the  final  stage  before  the  application  to  the  multi- 
target  tracking  problem.  The  purpose  is  to  get  an  actually  implementable 
algorithm,  that  connects  the  results  of  the  last  chapter  on  the  decompo¬ 
sition  of  the  estimation  problem,  to  the  results  of  Chapter  7  on  the  es¬ 
timation  of  a  single  Gauss-Markov  process. 

10.1  Problem  Statement 

The  system  considered  here  is  the  one  of  Chapter  2  Section  3,  under 
the  assumptions  (A.1]-(A.2}-(A.3)  of  the  last  chapter,  plus  an  additional 
Linear-Gaussian  assumption. 

As  in  Chapter  7,  an  algorithm  foT  estimating  the  discrete  or  hybrid 
state  sequence  of  the  multiple  system  will  use  some  pruning  rules  based 
on  the  distributions  P  (jc (k) , dk|  and  ? (x 00. , dk i Yk) ,  where  now  x(k)  « 

{x100,...,XpQ0)  and  dk  ■  {dk,...,dk}.  However,  under  the  Linear- 
Gaussian  assumption,  the  expressions  foT  these  distributions  arenot  as 
simple  as  for  the  single  Gauss-Markov  process  of  Chapter  7.  We  have: 

P0cCk),dk|Yk)  •  PCd^Y*)  Z  PCxCk)[dk,Yk,Ak)  x  P(Akidk)  (J0.1) 

(A1*) 

P(x(k) | dk,Yk,Ak)  is  the  conditional  density  for  the  multiple  state  xCk) 

k  k 

on  a  discrete  state  trajectory  d  and  a  data  hypothesis  sequence  A  up 

V 

to  time  k.  The  dynamics’  independence  implies  that  on  a  given  d  and 
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P(x(k)idk,Yk,Ak)  -  n  P(x.  (k) 1 dk,Yk,Ak)  (10.2) 

i-1 

Each  distribution  P(x^(k)jd  ,Y  ,A  )  is  a  Gaussian  density,  whose  mean 
*k(k^dk'Ak)  and  covariance  ^j]tCdk,Ak)  depend  on  both  dk  and  Ak,  and 
can  be  determined  recursively  from  the  Kalman  Filtering  Equations  (7.2)- 
(7.7),  where,  at  each  time  k,  the  state  estimate  is  updated  with  the 
measurement,  if  any,  given  by  A(k) . 

Therefore,  P(x(k)| dk,Yk,Ak)  can  be  written  as  N[x(k)  -  k(dk,Ak) , 
Pk|k(dk,Ak)J  where: 


%|kCd1C.Akl  • 


ad  Pk,kCdl  A 


*J|k  «“•**> 

L^ikAk> 

o  '  Pj^CdVlJ 

and  P(x(k),dk|Yk)  -  P(dklYk)  Z  P(Ak!  dk)  *  N[x(k)  -  k(dk,Ak) . 

{Ak} 


(10.3) 


Similarly,  we  have: 

P(x(k),dkfYk)  -  P(dk[ Yk)  Z 

(Ak> 

N(x(lc)  -  Xkjk(dk,Ak),P 


/C2ff)hk|Pk!k(dk,Ak)| 

PCAkUk)\CJ,)pnklpk|k(dkiAk)[,< 

k|kCdk,Ak)]  (10.1) 


klk  k  k  k 

where  P  1  (d  ,A  )  is  the  covariance  of  the  state  sequence  x 

k  k 
on  d  ,A*. 


conditioned 


So,  for  both  algorithms  the  pruning  rules  in  the  gTaph  would  be 
based  on  a  functional  inequality  between  two  linear  sums  of  Gaussian 
distributions.  Unfortunately,  no  simple  test  can  be  found  that  is  equi¬ 
valent  to  such  an  inequality.  However,  some  priming  rules  based  on 
PCxOO,dk,Ak!Yk)  or  P0t0c),dk,Ak|Y^  would  be  easily  implementahle  using 
the  Lemma  of  Chapter  7.  Thus,  instead  of  estimating  the  hybrid  or  dis¬ 
crete  state  sequence,  we  can  think  of  estimating  the  pair  (d^A^)  or 
1C  jc  fC 

Cs  ,A  j .  The  resulting  d  or  §  are  taken  as  our  state  sequence  estimates. 

10.2  Suboptimal  Algorithm  for  Estimating  the  State  of  a  Multiple 
Gauss-Markov  Process 

In  this  section,  we  develop  an  algorithm  for  estimating  the  pair 
(d'.A’S  or  (s^A^)  of  our  multiple  system.  The  derivation  is  very  close 
to  the  one  of  Algorithm  1  or  2  of  Chapter  7. 

We  define  the  pair  CdOO.AQO)  as  being  a  hypothesis  at  time  k  on 
the  multiple  system.  A  hypothesis  tree  can  be  associated  with  our  system. 
In  this  tree,  each  node  corresponds  to  a  distinct  hypothesis  h(k)  at 
time  k  and  each  branch  represents  a  transition  to  some  new  hypothesis 
at  the  next  instant  of  time.  (See  Fig.  10.1). 

k  k 

According  to  the  last  section,  the  distributions  P(x(k) >h  | Y) 
and  PCxCk),hk|Y]t)  are  equal  to: 


PCx(k).hkiYk) 

P(x(k),hk|Yk) 


PChk|Yk)  x  N[*(k)  -  ^c[kChk)>Pk|kchk)3>  CIO -5) 


P(hk|Yk)  x 


t/c^)nPlPk[kChk)l 

>  C2»)npk| Pk 


(27T)npk|Pklk(hS  | 


xN[,(M  -%!kCHh,Pklk(hk)] 


(10.6) 
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The  a  posteriori  probability  P  (hk|  Yk)  is  equal  to  PCh^.Y^/PCY^] 
where  P(Yk)  is  a  normalization  term;  and  P(hk,Yk)  can  be  determined  re¬ 
cursively  from: 

PCh*,^)  *  P (A GO |d GO]  *  PCda)|dCk-l))  x  PCYGO|h.k,Yk'1)  x 
P(hk"1,Yk"1]  (10.7) 

P(A(k)|d(k))  and  P(d(k) | d(k-l))  are  given  by  the  model  assumptions. 

k  k-1  k 

P(.Y(TO|h  ,Y  )  is  the  innovation  density  on  h  for  the  multiple  system. 

We  have,  using  Chapter  7: 

P(Y(10  |hk,Yk_1)  >  l  P(Z- Oc)  (dk,hk"1,Yk_1)  x  P(Y  CA(k)) 
i*l 

P  « 

*  P(YfCACk))  x  N^OO  -  Gi(di(k)  ,k)x^j ^(h  ) , 

B.(hk)J  (10.8) 

where  B.(hk)  «  ^(d.Gt)  ,k)pJjk_1Chk)H^(di(k)  ,k)  ♦  R.(d.Ck),k)  (10.9) 

K  K  K 

An  algorithm  for  estimating  the  multiple  sequence  h  or  (h  ,x  ) 
would  use,  as  for  Algorithm  1  or  2  of  Chapter  7,  some  pruning  rules  in 
the  hypothesis  tree  based  on  P(x(k) ,hk[Yk)  r  P(x(k) ,hk[ Yk)  which  are 
Gaussian  distributions.  To  get  a  simple  implementation  of  these  pruning 
rules,  we  can  now  apply  the  Lemma  of  Chapter  7. 

In  the  following,  we  only  derive  the  pruning  rules  and  the  cor¬ 
responding  algorithm  for  estimating  the  discrete  state  sequence  of  the 
multiple  system.  The  derivation  for  the  hybrid  state  estimation  algo¬ 
rithm  is  simply  obtained  by  formally  replacing  P^j^O^)  by  Pk^k(hk) 
and  g  by  y. 
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We  get  immediately  the  following  pruning  rule  MRl(.L),  that  compares 

If  V 

two  sequences  and  h 2  terminating  at  the  same  discrete  state  d(k)  at 
time  k  in  the  hypothesis  ttee: 


V  V 

Pruning  Rule  MR1(1Q:  Let  h^  and  h^  be  two  hypothesis  sequences  such 
that  d^(k)  *  d^Ck) •  Let 


k, ,  -1 


8Ch*,h!p  -  ^  o£|lca$ 

i  k  i  k  PO'jIV1'] 

«i|k°9  -^ik01!”  *2tn^77  * 


Z  in  -  - 

l-i  iPtifcOpI 


where  *Pk)  *Pk|  3X9  the  estimates 

v  k 

and  covariances  of  subsystem  i  on  h^  and  h2<  If,  for  all  i*l,2,...,p, 
^jk^l^  <  Pk|k^2^  31,41  S(h^,h^}  <_  0,  then  h^  can  be  immediately  dis¬ 
carded  at  time  k,  without  increasing  the  probability  of  an  error  in 
estimating  the  sequence  h  . 


Proof:  The  proof  comes  directly  from  the  Lemma,  the  expression  for 
PfrOO.h^Y*)  and  the  fact  that  PCACk)  [  d(k))  depends  only  on  the  discrete 
state  at  time  k. 

As  in  Chapter  7,  there  exists  a  recursive  relation  on  the  coeffi¬ 
cient  3,  that  is  independent  of  the  measurements.  This  relation  leads 
also  to  another  pruning  rule  as  follows: 
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Pruning  Rule  AMR1(L):  Let  h*  be  an  hypothesis  sequence  up  to  time  k.  If 

K  k 

for  all  hQ(k+l) , . . . ,hQCK) ,  there  exists  a  sequence  h‘  *  (h  ,hQ (k+1) , . . . , 
hg(K)),  where  the  corresponding  dCk)  is  equal  to  dg(k),  such  that: 


k  k  ?  K  (Bi^i 

1  E  E  in 

0  i-1  jak+l  !BiChJ)l 

p  |pit  t0*^! 
i=l  |P*.*Gn| 


PJUCHk)| 


♦  r  in  r*  ^ 

i-l  |Pk)k0^ 


then  h^  can  be  immediately  deleted  at  time  k,  without  increasing  the 
probability  of  an  error  in  estimating  . 

Using  MR1(L)  and  AMR1C.L),  an  algorithm,  similar  to  Algorithm  1  of 
Chapter  7,  can  be  derived  for  estimating  the  sequence  Cd^.A^]  of  the 
multiple  system. 

However,  we  must  also  use  the  decomposition  possiblities,  introduced 
in  the  last  chapter,  that  can  lead  to  a  substantial  reduction  in  the 
amount  of  confutation.  If,  for  all  k  ■  0,1,..., K,  the  sets  Yfk)  R^(k) 
defined  in  the  last  chapter  are  mutually  exclusive,  then  we  can  solve 
the  estimation  problem  by  running  p  independent  algorithms.  For  each 
algorithm  i,  a  hypothesis  tree  for  the  corresponding  subsystem  i  is 
generated.  Then,  some  pruning  rules  based  on  PCx.j.Ck) »h^| Y^)  are  used 
to  delete  at  time  k  some  hypothesis  sequences  h.  during  the  formation  of 
txee  i.  Using  the  same  arguments  as  in  Chapter  7,  these  pruning  rules 
can  be  stated  as  follows: 


k  k 

Pruning  Rule  R1(L)  in  tree  i:  Let  h  .  and  h  .  be  two  hypothesis  se- 

l » i  2 ,  i 

quences  up  to  time  k  in  tree  i  such  that  the  corresponding  d.  . (k)  is 

1)1. 

equal  to  d.  . (k) .  Let 

2,1 
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Pk|k<h^‘‘«i|k<i>  -  ^Iktbj.i”  *  2ln^H‘Y  ’ 


Pth^iin 


pklk<1)l 


lPklk®2,i^ 


If  Pk|kth![,i)  <  pkik^.t>  ™i  th“  hi.i  “> be  i-edi- 

ately  discarded  at  time  k,  without  increasing  the  probability  of  an  error 

in  estimating  (d^,A^  in  tree  i. 


in  tree  1. 


We  have  also  the  additional  pruning  rule: 


Pruning  Rule  AR1CI0  in  tree  i;  Let  hQ  ^  be  a  hypothesis  sequence  in 

tree  i  up  to  time  k.  If,  for  all  hn  .  Ck+1) , . . • >h  . (K) ,  in  tree  i, 

K  k 

there  exists  »  Chi  ,kQ  iCk+l), . . .  ,hQ  ^K)),  where  the  corresponding 
d.(k)  is  equal  to  dQ  iQc),  suchj&jat  P^C^q^)  <  P^jk^)  and 

k  K  ,31  iPkuChJll 

0Ch*  ,h.j<  E  An  -  1-  ♦  in  -  | - 

°'1  1  ;i-k,i  i.,®*)!  ip^o^i 

!pK|K0>f>l 

][ 

then  hn  .  can  be  immediately  discarded  at  time  k,  without  increasing  the 

0,1 

K  k 

probability  of  an  error  in  estimating 

If  the  independence  among  the  algorithms  can  be  preserved  until 
time  K,  then  the  estimate  h  is  simply  given  by  (h^, . . .  ,h")  resulting 

*  r 

from  the  subalgorithas. 

If,  at  a  given  time,  some  overlaps  occur,  then  we  have  to  go  back 


to  a  combined  algorithm  using  MR1(L),  AMR1(L). 
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K  K 

The  algorithm  MA1(L},  for  estimating  the  pair  (d  ,A  )  of  the  mul¬ 
tiple  system,  uses  the  combination  of  the  local  algorithms,  whenever 

k  k 

it  is  possible,  and  the  combined  algorithm.  We  denote  byj*^  and  Jf' 
the  remaining  hypothesis  sequences  at  time  k  in  tree  i  and  in  the  global 
tree  respectively.  The  algorithm  can  be  stated  formally  as  follows: 


ALGORITHM  MAI  CM 


STEP  0: 


*Si.i<y  ■  4  ■  poi-itv  ■  p5'  -  R°- for  a11 

i*  1*2, ... ,p 

Test:  V  Ci , j  1  6  U,p]2,  YCO)  A  R°  A  R°  =  0 


YES:  Go  to  SUBALGORITHM  SA1 (L) 
NO:  CONTINUE 


STEP  k:  ( 

I  POtk’1|Yk-1]*4.1|k.1Chk-1h  P^ilk.!^'1) 

for  all  i  *  1,2, ... ,p, 
for  all  hk“X  6 
RECEIVE  Y(k) 

COMPUTE  for  all  hk  6  x  all  the  possible  multiple 

hypotheses  at  time  k,  P(h.k|Yk),  Pk[kfrk)  • 

APPLY  MR1(.L) 


APPLY  AMR1CL) 

STORE  the  remaining  sequences  in 
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k  »  k+1  and  repeat  until  k  »  K 

STEP  K:  | 

|  P(hK|YiC)  for  all  hK  6  J<rK  . 

Find  Sup  PCh^jY*)  .  The  corresponding  sequence  is 

hKe^ 

equal  to  h^.  STOP. 


SUBALGORITHM  SA1 CL) : 

STEP  k:  for  all  i  «  1, 


PCh^1^"1 

pk-iU-ithJ'1)-  for  a11  "i'1 6 


*k-l 


k-l^i  5 


l 

COMPUTE  for  all  i  *  l,...,p,  R^k) 

(A  method  for  computing  R-^k)  will  be  derived  for  the 
multitarget  application) 

TEST:  for  all  (i,j)  6  {l,p]2,  Y(k)  n  R^k)  n  (k)  *  0 
YES :  Continue 

NO:  Go  to  Step  k  of  Algorithm  MAI  CL). 


RECEIVE  YCk) 

k  k-1 

COMPUTE  fdr  all  i  ■  1,2,. ..,p,  for  all  h.  *  all  the 

i  k  i  k 

possible  hypotheses  in  tree  i  at  time  k,  S^Gu),  P]c[  * 
P(h*|Y*) 


APPLY  R1CL)  in  tree  i,  for  all  i  *  l,...,p 
APPLY  AR1CL)  in  tree  i,  for  all  i  »  l,...,p 


STORE  the  remaining  sequences  in 
k  ■  k+1  and  repeat  until  k  »  K 


ye*. 
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STEP  K:  For  all  i 


i  *  l,...,p(.#f 


PJ-x 

(  PCh^Y^  for  all  hK  6 


r  l 


Find  Sup  P(h^|Y^)*  The  corresponding  sequence  is 
1  1 


our  estimate  h.  . 
Then  fiK  »  (h^ , . . . ,hp) 


END. 


So,  Algorithm  MAI  CL)  associates  the  pruning  rules  coming  from 
Chapter  6  and  7  with  the  decomposition  result  of  Chapter  9.  A  compu¬ 
tation  of  the  measurement  regions  R^0O  will  be  done  in  the  next  chapter 
for  the  application  to  multitarget  tracking. 

A  similar  algorithm  MA2C.L)  can  be  derived  for  the  estimation  of 


av  AJ^  jr 

The  resulting  estimates  d  or  Cd  ,x J  given  by  the  most  likely 
(d^A1^  or  (dK,AK,xK)  are  not  necessarily  the  MAP  estimates  of  dK  and 
(d^x*)  respectively. 

The  continuous  state  estimates  are  computed  for  the  most 
likely  data  hypothesis,  i.e.,  the  data  association,  and  do  not  take  into 
account  all  possible  data  hypotheses  weighted  by  their  probability. 

Reid  [9]  and  KeveTian  [10]  take  the  same  approach  in  their  multitarget 
tracking  algorithm.  So,  in  strict  estimation  terms,  our  estimates  for 
the  discrete  or  hybrid  state  sequences  are  optimal  only  in  a  rather 
unique  sense.  On  the  other  hand,  this  different  estimation  problem 
leads  to  much  more  tractable  computations.  Therefore,  we  suggest  this 
suboptimal  scheme. 

Using  Algorithm  MAI  CL)  or  MA2(.L),  the  storage  requirements  can 
still  be  too  large  with  respect  to  our  resources.  In  this  case,  the 
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algorithms  are  truncated  to  some  manageable  length  in  order  to  keep 
the  number  of  hypotheses  at  a  reasonable  level.  We  decide  upon  the  pair 

/V  A 

CdCk-K.^),  ACk^Kr])  at  time  k  by  choosing  the  most  likely  hypothesie  se¬ 
quence  from  k  -  K  -  1  to  k.  The  decomposition  property  applies,  when¬ 
ever  the  sets  Y0c)  A  fLQO  are  mutually  exclusive  during  an  interval 

[k.  -  Kr,  k],  even  though  there  has  been  some  overlap  before  time  k  -  Kr> 

The  proposed  algorithm  in  the  next  chapter  for  solving  the  multi- 

target  tracking  problem  is  based  on  a  MA.1CL)  type  algorithm  truncated 
to  a  length 

10.3  Summary 

In  this  chapter,  an  algorithm  for  estimating  the  discrete  or  hybrid 
state  sequence  of  a  multiple  Gauss-Markov  process  has  been  presented. 

It  contains  both  the  implementable  pruning  techniques  of  Chapter  7  and 
the  decomposition  possibilities  of  Chapter  9. 

In  the  next  chapter,  we  discuss  an  example  application  to  multi- 
target  tracking. 
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CHAPTER  11 

APPLICATION  TO  MULTITARGET  TRACKING 


The  purpose  of  this  chapter  is  to  apply  the  algorithm  presented  in 
the  last  chapter  to  the  multitarget  case.  In  particular,  a  computation 
of  the  measurement  regions  R^QO  will  be  done,  and  the  decomposition 
possibilities  will  be  discussed  in  this  context  through  a  simple  numeri¬ 
cal  example. 


11.1  Model  and  Objective  for  the  Multitarget  Tracking  Problem 


The  model  for  the  set  of  targets  has  been  presented  in  Chapter  2  as 
an  example  of  a  multiple  Gauss-Markov  process.  We  recall  briefly  the 
main  features. 

The  system  considered  consists  of  p  targets.  We  assume  that  each 
target  goes  successively  from  an  unborn  discrete  state  U,  to  an  alive 
state  A,  to  a  dead  state  D.  Planar,  straight-line  motions  are  assumed 
for  the  targets  (see  Fig.  11.1).  So,  the  state  s^(k)  *  { d . (k)  ^ GO } 
is  defined  by: 


d.(k)  6 


x.Ck) 


1U,A,D}_ 

xjoo" 
vjoo 
xjjoo 

v?(k) 


if  d.00  -  A 


Each  target  i  can  be  detected  by  the  surveillance  system  with  a  con- 
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stant  probability  of  detection  P^.  The  probability  for  a  false  alarm 
is  a  constant  P^.  We  assume  that  the  detectability  of  each  target  is 
independent  of  the  otheTS  and  that  the  process  generating  a  false  alarm 
is  independent  of  the  other  false  alarms.  We  assume  also  independent 
target  dynamics.  So,  we  have  for  each  target  i: 

-  the  discrete  dynamics: 


the  continuous  dynamics:  if  d^Qc)  *  d.Qoll  *  A 


xiCk+l] 


1  1 
0  1 
0  0 
0  0 


0 

0 

1 

0 


*iW  +  vtG0  ,  *^00  ~  NCO.cy 


xiC0)  -  NC0,P“) 


the  measurement  equation:  if  d^OO  *  A, 


tjOO 


XjOO  ♦  viCh)  ,  -  NCP.R.); 


10  0  0 
0  0  10. 

if  d-^k)  6  {U,D},  no  measurement  is  generated;  so,  we  can  set 
ti(k)  to  0. 


-  the  false  alarms  distribution: 


2.  (k)  -  N(0,Z)  ,  Z  >  0 

f 


-  the  data  hypothesis  statistics: 

m  j(k) 

PtACk)|d00)  »  P/  x  Cl-Pd) 


M  (k)  0  MfCk)  , 

*"v  ‘  P,  '  ,  where: 


’m' 


x  *f 
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M^OO  *  *  detections  under  A(k) 

M^OO  *  #  missed  detections  under  A 00 
M^OO  ■  #  false  alarms  under  AQ0 

The  purpose  of  the  surveillance  system  is  to  decide  upon  the  dis¬ 
crete  states  of  the  targets,  and  to  estimate  their  locations  and  velo¬ 
cities,  with  a  maximum  delay  Kr,  while  minimizing  the  probability  of 
an  error  in  the  decisions.  The  delay  Ky  takes  into  account  both  the 
computational  resources  and  the  necessity  for  a  timely  decision. 

Algorithm  MAI  CL}  of  the  last  chapter,  truncated  to  a  length  K^, 
applies  directly  to  this  problem,  since  it  computes,  at  each  time  k  and 
for  each  target  i,  the  discrete  state  d.(k  -  ,  the  data  association 

A 

A^k  -  Kr),  The  location  and  velocity  estimates  can  be  simply  given 

i  A 

by  the  continuous  state  estimate  i  A^(k-Kr))  which 

r1  r 

is  stored  at  time  k  -  Ky  in  the  algorithm.  Using  MAI  CL) ,  we  minimize 
the  probability  of  an  error  in  deciding  upon  the  discrete  states  of 
the  target  and  in  performing  the  data  association. 

To  take  advantage  of  the  decomposition  possibilities  in  this  al¬ 
gorithm,  we  must  compute  explicitly  the  measurement  regions  R. Ck} • 

This  would  show  the  optimality  of  a  gating. 

11.2  Optimality  of  a  Gating 

We  consider  a  given  target  among  the  set  of  targets,  whose  exis¬ 
tence  is  assumed  on  at  least  one  remaining  hypothesis  in  the  tree  at 
time  k-1.  For  each  such  hypothesis,  we  have  an  estimate  of  the  position 
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and  velocity  of  the  target  and  the  corresponding  covariance  (see  Fig. 
11.2).  So,  the  predicted  values  for  the  next  measurements  coming  from 
that  target  and  their  covariances  can  be  computed.  If  these  covariances 
are  small  and  we  receive  at  the  next  instant  of  time  a  measurement  far 
away  from  these  predicted  values,  it  is  unlikely  that  it  should  be  asso¬ 
ciated  with  this  target,  since  it  would  require  a  very  unlikely  value  of 
the  measurement  noise.  Rather,  we  would  suspect  that  this  measurement 
is  a  false  alarm  or  is  due  to  another  target.  This  intuitive  approach 
is  taken  in  £9] -{10],  where  an  arbitrary  gate  is  drawn  around  the  pre¬ 
dicted  value  of  the  measurement  for  each  hypothesis;  a  measurement 
falling  outside  this  gate  is  not  associated  with  the  target.  A  justi¬ 
fication  for  this  gating  can  be  found  in  our  optimal  hypothesis  deletion 
techniques  as  follows. 

We  have  to  look  for  a  hypothesis  at  time  k,  that  consider  the  re¬ 
ceived  measurement  z(Tc)  as  a  false  alarm,  and  could  delete,  using 
R1  (L) ,  the  measurement  association  of  z(k)  to  the  target,  for  some  z(k). 
These  values  of  z(k)  define  a  region  outside  which  no  measurement  asso¬ 
ciation  is  to  be  made  to  the  target.  We  can  think  of  simply  extending 
each  hypothesis  sequence  on  the  target  at  time  k-1  by  the  two  following 
hypotheses:  the  target  is  still  alive  and  the  measurement  z(k)  is  as¬ 
signed  to  it,  or  the  target  is  alive  but  z(k)  is  a  false  alarm  (see 
Fig.  11.3).  Unfortunately,  it  can  be  shown  that  there  exists  no  measure¬ 
ment  for  which  the  first  extended  hypothesis  can  be  deleted  by  the 
second  extended  hypothesis,  using  Rl(.l):  the  corresponding  comparison 
coefficient  3  is  equal  to: 

p 

iT00£'1i(M  •  2  *  *«  fj|f  • 


(a,  a. 

a&um.) 


>  (A,  i(4.)  jut)  aut^/tuid. 
to  fongd-t ) 

(A,z(&)m>  x  faJoe 

aSattm. ) 


*  (Ai  M*) 

(A,  z(4)ua  x  S&e 

aJowr*) 


Jo  oj» (,^/t^d- 
h  hxr^el*) 
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which  is  positive,  for  all  z(k). 

Instead,  we  must  try  to  compare  each  hypothesis  at  time  k  that 
associates  a  measurement  z(k)  with  the  target  against  another  hypothesis: 
that  a  secondary  target  is  just  bom  at  time  k  hut  not  detected  and  zCk) 
is  a  false  alarm.  For  this  alternative  hypothesis,  we  assume  only  a 
priori  information  on  the  secondary  target  state;  an  a  priori  mean  and 
an  a  priori  (large)  covariance,  i.e.,  no  information  at  all  (see  Fig. 
11.4).  This  hypothesis  implies  that  the  preceeding  measurements  asso¬ 
ciated  with  the  primary  target  are  now  considered  as  false  alarms.  If 
we  assume  that  once  in  the  past  the  existence  of  the  primary  target 
has  been  proved,  because  all  the  survivors  went  to  the  alive  state  A, 
then  the  alternative  hypothesis  is  only  considered  for  our  gating;  we 
do  not  have  to  take  into  account  the  hypothesis  that  the  secondary 
target  is  just  bom  and  that  a  measurement  is  associated  to  it.  This 
argument  applies  to  the  single  object  case  as  well. 

It  is  possible  to  determine  the  set  of  measurements  z(k)  for  which 
all  the  hypotheses  that  would  associate  the  observation  can  be  deleted 
by  that  last  alternative  hypothesis.  This  corresponds  to  our  measure¬ 
ment  region  R^(k)  for  target  i  at  time  k.  An  actual  computation  of  such 
a  region  proves  the  existence  of  an  optimal  gating.  It  can  be  done  in 
the  following  way. 

Let  us  consider  a  specific  target,  the  primary  target,  at  time  k-1, 

^-1  k-1 

and  the  corresponding  remaining  hypothesis  sequences  Jtr  =  Ch^  , 

1  <_  j  <_  r)  (the  index  i  of  the  target  is  omitted  for  simplicity  in  the 

k-1  Jk-l 

following  notations)  terminating  at  d(k-l)  *  A.  For  all  h.  6  Jc  , 

k-1  k-1 

we  have  an  estimate  1  a  covariance  pjt_^  J  k-1  ^ ' 


,(A,  z(4)io  Moi&nai  if  fie 

fotyaf-) 


^  (A,  t(A)ua  <t  JUoe  ahxm) 
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We  can  compute  the  predicted  value  z ^  (k)  of  the  measurement ,  corresponding 
to  h^-1,  from  2^00  a  H\jk.-l^j^’  ^  t*ie  covariance  B(hj)  =  HPkjk_1(h^HT 
+  R,  for  all  h^  *  Chj'^.h^OO)  where  h^Ck)  is  the  hypothesis 
CdQtl  ■  A,  the  measurement  originates  from  the  primary  target) .  Under 
our  model  assumptions,  we  have: 


POyY*)  •  POv^.Y*-1)  x  Cl  -  P1(J)  *  Pd*  N[z  (k)  -  z.  (k) , 


Cll.l) 


^kjk^j^  are  uP^ate<^  with  zCk)  using  the  Kalman  Filtering 
Equations  (7 . 2] - (7 . 7) . 

k. 

We  compare  each  hypothesis  h.  with  the  alternative  hypothesis 

]c 

hQ  »  (the  secondary  target  is  just  born  at  time  k,  not  detected  and 
z(k)  is  a  false  alarm).  We  have: 


PChg.Y11*1)  =  PCh^.Y^1)  x  PQ1  x  Cl  -  Pd)  ^  Pf  x  N{z(k),rj 

CU*2) 


k  k 

*k|k^(P  *  xo>  Pk|k^V  "  P0  is  the  a  Priori  "larSe"  covariance, 
v  v 

We  have  P^^Ch.^)  <  Pq,  since  under  h^  we  have  updated  the  estimate 

and  covariance  at  least  once.  So,  according  to  R1(L),  we  don’t  assign 

k  k 

z(k)  to  the  target  if  the  corresponding  S(hg,h.)  is  less  or  equal  than 

k-i 

zero.  We  denote  by  R(Tk  )  the  set  of  measurements  of  the  space  Z  for 
which  S(hQ,h^)  >  0,  i.e., 

RQk'1)  *  {z(k)  £  Z  s.t.  S(hg,h^)  >  0} 


Then,  the  measurement  region  RCk)  corresponding  to  the  target  is  simply 
given  by: 
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R(k)  *  U  R-Qv^'1) 

*JS  J 

In  Appendix  D,  we  show  that, 


zGO  S  ROk-1)  <  *•>  SCh^hh  >  0 

•  j  u  J 


where 


<->  CzGO  -  m.Ck))lE.CRKz0c)-m.0c))  <  EjO) 


(11.3), 


EjCH  ■  Qt  ♦  HlP^nO^  -  P^'V]-1  -  S'1  Cl  1.4) 

mjOO  ■  E^CWCR  +  HIp;lk-iOiJ}  -  Cll.S] 

=)«  -\|k-l«£  *  pk|k-l4rlt5lk|ic-i^  - 

x0)  C11.63 

e.Ck)  =  CHc.(k)]Tlz  -  R  -  HCF^^Oi^  -  P^W]’1  CH^ 00) 


♦  5  ..00 

p.  Cl-P10) 

+  ^  *nfl-p  )P  *  ^  p 
u  dJ  f  01 


'■r\ 

+  in  (if* 


(11.7) 


«jW  ■  C»k|k.i»j)  -  xo)T"’o  -  \|k-la‘J»1®k|k-lthj>  -  V 


PQk"1  Y^"1) 

♦  2  in  - 1—. - e-T—  JtTl  — 

PChg  ,Y  )  |P 


’kjk-l^j^ 


C11.8) 


We  can  assume  that:  P^  >  (.l-P^P^,  E  >>  R»  Pgj  +  P^g  <  1*  $j(k) 

It  k,  k  k-1 

is  positive,  since  it  corresponds  to  S(h  ,hg)  where  *(h^  ,h.^(k)) 

k  k  1  v 

Cif  3(h*,h^  j  <_  0,  then  we  could  immediately  discard  fu ,  whatever  t(k) 

is).  So,  if  H[P“Jk_1(hJ)  -  P ~ 1 ] ' 1HT  <  I  -  R  (11.9),  then  E . (k)  >  0  and 
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t  (k)  >  0.  Condition  01.9)  is  satisfied,  if  Z  is  much  larger  than  R  and 
k 

^kjk-l^Lj^  nuc^  smaller  than  P^,  which,  is  generally  the  case. 

k.-l 

In  this  case,  R(h^  )  defines  an  ellipse  centered  at  m^Ckl-  The 

measurement  region  R0t)  fo r  this  target  is  the  union  of  several  ellipses 
k-1 

R(hj  )/  whose  computations  are  quite  simply  using  0.1  -33 -01 •  81.  An 

approximation  of  R(k)  without  loss  of  estimation  performance  is  given 

by  one  ellipse  circumscribing  all  others  Csee  Fig  .  11.51.  Thus,  all 

the  measurements  falling  outside  this  region  R[k)  won't  be  associated 

with  the  target.  This  shows  the  existence  of  an  optimal  gating! 

A  special  case  of  interest  is,  when  the  distribution  of  the  false 

alarms  tends  to  a  uniform  distribution,  i.e.,  Z  -*•  ®  ,  and  P^  is  a  very 

large  initial  covariance,  so  that  we  can  assume  Pj,  ^  «  PQ  for 

k  k-1 

all  h ..  fi  *  x  hd00  •  In  this  case,  the  parameters  m^  QO  ,  E.  (k) , 
of  the  ellipse  take  the  following  simple  form: 


m.Ck)  t 

E.(k)  t  (R  ♦  HPk|k.1Ch^HT]'1  »  [BQVjlj”1 

Pj  l"^n  |r 

e.Ck)  s  2  in  *  2  in  ^  ♦  in  || 


01.10) 


01.11) 


*  5  j  CM 


(H*k|k.1(h'))Tr-1(Hiek|l.1(h^) 


01-12) 


K*  1 

The  ellipse  RCh^  )  becomes  centered  at  the  predicted  value  of  the 
measurement  (k) . 

The  previous  algorithms  of  Reid  [9]  and  Keverian  [10]  use  a  gating 
around  the  predicted  value  of  the  form:  (z(k)  -  i. (k))TCB(h^))'1(2(k)  - 
^(k))  <  n.  where  n  is  an  arbitrary  threshold.  Using  the  above  result, 
we  now  have  a  mean  of  determining  the  optimal  threshold  equal  to  e.fk) 
defined  by  01.12) . 


So,  for  each  target  i  at  tine  k-1,  we  compute  the  ellipses  R^(k) 
using  .  If,  for  each  k,  the  sets  Y(k)  D  R^Qc)  are  mutually 

exclusive,  then  our  multitarget  problem  is  equivalent  to  p  independent 
single  target  problems. 

We  now  look  at  a  simple  numerical  example  to  illustrate  the  appli¬ 
cability  of  this  optimal  gating. 

11.3  Example 

The  following  example  focuses  on  the  decomposition  possibilities 
in  the  multitarget  tracking  algorithm,  i.e.,  computation  and  application 
of  the  optimal  gating. 

For  this  purpose,  we  consider  the  multiple  system  as  consisting  of 
two  targets.  We  assume  that  the  existence  of  these  two  targets  has  been 
decided,  because  the  single  remaining  survivor,  ot  the  most  likely  se¬ 
quence  in  case  of  a  truncation,  went  through  the  alive  state  for  both 
targets  at  one  time.  We  assume  also  that  the  covariances  for  the  con¬ 
tinuous  state  estimates  have  reached  their  steady-state  values  Cin  5 
to  10  steps  with  the  following  values).  For  simplicity,  at  each  time 

k  we  compute  for  each  target  only  the  ellipse  corresponding  to  the  com¬ 
ic. 

parison  of  the  hypothesis  sequence  h  . C  primary  target  i  always  de- 

A»i 

tected  and  a  measurement  is  assigned  to  it)  with  the  hypothesis  sequence 

h^  ^  Csecondary  target  i  is  just  born  at  time  k,  but  not  detected  and 

the  measurement  is  a  false  alarm).  We  study  the  variation  in  the  sice 

of  this  ellipse  as  a  function  of  the  process  and  measurement  noises 

and  the  efficiency  of  the  gating  as  these  two  targets  are  crossing. 

We  consider  Q,  and  R.  matrices  of  the  form:  for  i  *  1,2, 
xi 


0  0  0  0 


0  q  0  0 
0  0  0  0 
0  0  0  q 


r  Q 


0  r 


where  the  process  noise  enters  only  the  velocity  components. 
We  take  the  following  values  for  the  other  parameters: 


-  discrete  dynamics:  =  P  ^  »  0.1 


-  initial  conditions: 


100  0  0  0 

0  5  0  0 

0  0  100  0 

0  0  0  5 


*  detection  and  false  alarm  probabilities: 

P,  *  0.66;  Pr  *  0.33 
d  f 

-  false  alarms  covariance: 


0  500 


for  i  * 


Assume  that  at  time  k-1,  we  have  the  following  probabilities  and 
estimates  for  target  1  and  2: 


pCh^'!|Yk“1]  -  PCh^l^'1)  .  O.S 

PCho’ll^”1)  *  P^O^^  *  0,01 


1  ..k-1  I  *5  2  ..k-1.  a  -5 

*k|k-l  1  ,n  ^kj  k-1  ^A,2 
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We  compute  the  ellipses  Cin  fact  circles)  for  different  values  of  q  and 

i  k 

r.  Since  we  take  r«  500  and  we  have  P^j  «  P^,  these  circles 

are  centered  at  the  predicted  values  ("^q)  ^0T  target  1  and  ^or 

target  2.  In  Fig.  11.6  we  have  plotted  the  radius  p  of  these  circles 
as  a  function  of  q  and  r. 

We  see  that  for  q  or  r  greater  than  1,  the  radius  grows  rapidly, 

so  that  the  gating  becomes  useless.  As  r  tends  to  0  for  q  equal  to  1, 

this  radius  tends  slowly  to  ■*•,  because  of  the  term  Jin  j-|p  .  For  r=l , 

i  k 

as  q  tends  to  0,  it  tends  also  slowly  to  «*>  because  of  -  Jtn|  P^j  ^)  j 

in  5  j  00  . 

Thus,  we  see  that  the  size  of  the  gate  is  relatively  insensitive 
to  noise  modelling  errors  for  small  values  of  q  and  r  Cl ess  than  0.1  in 
our  example) . 

Now  for  q  ■  r  ■  1,  we  study  the  efficiency  of  the  gating  as  two 
targets  are  crossing.  To  limit  the  number  of  hypotheses,  and  subsequent¬ 
ly  the  number  of  circles,  we  assume  *  1. 

*1  k 

We  start  at  time  k-1  with  the  above  estimates  1^ 

2  k 

‘Tc|k-l^hA  2^  hypothesis  probabilities  for  target  1  and  2.  We  re¬ 
ceive  some  sets  of  measurements  Y(k}  ,Y(k>l)  ,YCk*2) , . . . ,  and  each  for 
each  time  k,  k+1,  k+2,...,  we  compute  the  optimal  gates  foT  target  1 
and  target  2  and  apply  the  decomposition  theorem,  whenever  this  is 
possible  (see  Fig.  11.7-11.9): 


step  k:  Y(k)  A  R  (*)  r\  R^Ck)  -  0 

Assign  ZjCk)  to  target  1,  z^CXl  to  target  2, 
z,(k),  z.(k)  to  the  set  of  false  alarms. 

3  ‘  4 


Cfigure  ll-7a) 
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step  k+1:  Y(k+1)A  R,  (k+1)  A  R„Ck+l)  *  0  (figure  ll-7b) 

~  i  ^ 

Assign  z^k+1}  to  target  1,  z2(k+l)  to  target  2, 
z30c+l)  to  the  set  of  false  alarms 

step  k+2:  Y(£+2)  A  R2(k*2)  f\  R^+21  t  0  (figure  ll-8a) 

Data  association  uncertainty  for  z^(k+2},z2(k+2) . 

Assign  z3Ck+2y,z4(k+2}  to  the  set  of  false  alarms. 

step  k+3:  (figure  ll~8b) 

Since  Kr  *  1,  MAP  decision  rule  assigns  z2(k+2)  to 

target  1  and  z^Qt+2}  to  target  2  at  time  k+2. 

Data  association  uncertainty  for  z  (k+31 ,z_(k+3) . 

1  ■+ 

(figure  ll-9a] 

MAP  decision  rule  assigns  z2Qc+3)  to  target  1  and 
ZjCk+3}  to  target  2. 

Computation  of  the  corresponding  R^+4),  R2(k.+4). 

Y(k+4)  n  Rj 0t+4)  A  R2CX+4)  =  0  Assign  z^k+4)  to 

target  1,  z2(k+4)  to  target  2,  z3(k+4) ,z4(k+4) ,zg(k+4)  to 
the  set  of  false  alarms. 

Y(k+5)  A  RjCk+5)  A  R2(k+S}  -  0  (figure  ll-9b) 

Assign  z2(k+5)  to  target  1,  z^k+5)  to  target  2, 
z3(k+S)  to  the  set  of  false  alarms. 

We  see  that  the  crossing  of  the  two  targets  introduces  a  correla¬ 
tion  at  time  k+2,  k+3,  k+4,  but,  after  that  time  the  decomposition 
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property  is  regained.  Therefore,  we  must  run  one  algorithm  for  both 
targets  only  for  steps  k+2,  k+3,  k>4. 

This  example  highlights  the  efficiency  of  our  optimal  gating  in  a 
very  simple  case,  where,  for  simplicity,  we  compute  only  the  gates  cor¬ 
responding  to  the  targets  being  detected  and  we  take  the  delay  equal 
to  one. 

Complete  simulations  would  have  to  be  run  to  show  the  real  capabi¬ 
lities  of  our  optimal  pruning  rules  in  Algorithm  MA1(L),  and  the  effi¬ 
ciency  of  our  gating  in  more  complex  situations,  where  we  have  many 
targets  in  the  area,  we  take  a  delay  greater  than  one  to  decide  upon  the 
discrete  state  of  the  targets  and  the  data  associations,  and  we  consider 
all  hypotheses  for  each  target. 

11.4  Summary 

In  this  chapter,  we  have  computed  an  optimal  gating  for  the  multi- 
target  tracking  problem.  We  have  shown  that,  under  some  approximations, 
this  gating  took  a  form  similar  to  the  one  in  the  heuristic  procedure 
of  Reid's  algorithm  [9].  Then  we  have  derived  a  simple  example,  that 
illustrates  its  sensitivity  with  respect  to  the  noise  parameters,  arid 
its  efficiency  as  two  targets  are  crossing.  These  encouraging  results 
are  only  the  first  stage  in  an  evaluation  of  the  performance  of  our 
multitarget  tracking  algorithm. 


F16URS  M.l 
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CHAPTER  12 

CONCLUSION  AND  FURTHER  WORK 


The  purpose  of  this  thesis  was  to  set  up  a  common  framework  foT 
different  hybrid  state  estimation  problems,  and  to  present  optimal  al¬ 
gorithms  for  estimating  the  state  sequence  of  single  and  multiple 
Gauss-Markov  processes.  We  derived  optimal  techniques  for  reducing 
the  amount  of  computation,  and,  in  the  multiple  system  case,  we  proved 
that  clustering  can  be  optimal  under  certain  conditions.  We  also  com¬ 
puted  an  optimal  gating,  which  decreases  the  computational  complexity. 

Simple  applications  to  a  failure  detection  problem  and  a  multi¬ 
target  tracking  problem  showed  encouraging  results.  Due  to  the  lack 
of  time,  we  were  not  able  to  perform  more  complete  simulations  in  order 
to  evaluate  the  real  performances  of  our  algorithms:  efficiency  of  the 
pruning  rules  and  the  gating,  sensitivity  to  the  various  parameters. 

Our  hope  is  that  this  thesis  will  provide  a  basis  for  rational  fif 
not  optimal)  procedures  applicable  to  diverse  hybrid  state  estimation 
problems. 

A  further  area  of  research  would  be  to  apply  this  approach  to  a 
distributed  system,  where  one  single  decision  maker  cannot  handle  com¬ 
pletely  the  estimation  of  the  hybrid  system,  so  that  this  task  must  be 
divided  among  several  local  decision  makers,  that  limit  the  communications 
between  them  at  the  minimum  level  required  to  insure  a  coordination  in 


their  decisions. 
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APPENDIX  A 

Some  Useful  Matrix  Identities 

Let  and  B.,  be  two  n  *  n  square  invertible  matrices.  We  have  the 
following  equalities: 

B"1  -  -  B^1)"^1  =  CBX  -  B^*1  CA.l) 

B"1  +  B^CB-1  -  B'V^j1  =  CB,  -  B^'1  (A. 2) 

In  Jazwinsky  [4],  we  have  the  following  equalities: 

T 

Let  P,R,M  be  n  *  n,  m  *  m,  m  *  n  matrices  where  P  *  P  >_  0  10(1 
T 

R  »  R  >0.  Then: 

Cl  ♦  PMTR_1M)_1  =  I  -  PMT(MPMT  ♦  R)_1M  (A. 3) 

(I  +  PmVHi)"1?  •  P  -  PMTCMPMT  ♦  R) _1MP  (A. 4) 

(I  +  PMTR_1M)“1PMTR  =  pmtcmpmt  +  R) 


(A. 5) 
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APPENDIX  B 


Proof  of  the  lemma  in  Chapter  7 


v  x  e  x,  fxcx)  i  f200 


CB.l) 


<->  Inf  2  (Jin  f,(x)  -  Inf.Cx))  >  0 
x€X  1 

<*>  Inf  [xT(B"1  -  B'^x  -  2(x  B11  ‘  x2B2  )x  +  X1B1  X1  “  X2B2  x2  + 
x€X  1  *  1 

°2 

2  sf  3  t  o 

The  extremum  is  equal  to:  xQ  *  (B^  -  B~1)”1(B”1x1  -  B^x^ 

It  is  a  minimum  if  B^  <  B^. 

Therefore: 


CB.l)  <=>  j  BL  <  B2 

j ■  oft1  -  xf21)(s2l  -  • 


l2x2>  *  2  ^  i  0 


<->  V  h 


2  .  Tr_-1  „-l,„-l  „-l, 


Tr  „-l 


2taSf+xitBi  -V (Bi “Vi  Bi  5xi  -  x2l-b2 


“z1'8!1  -  *  ^‘cb;1  -  BjhsJ1*!  >  0 


Using  (A. 1)- (A. 2) : 

(B.l)  <->  /2  in  ^  ♦  x^(B1  -  b2)"1x1  *  x]^  -  Bj)"^  _  2x2(B2  -  Bj)"^ 

V  B2 
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APPENDIX  C 

Proof  of  relations  (7. 15“)- (7. 16)  on  3  ,Y 


The  proof  is  straightforward  using  some  algebraic  manipulations. 

Let  d!^*1,  d!^*1  be  2  state  sequences  such  that  Ck)  *  d2(k)  and 
djCk+1)  =»  d2(k>l).  Let  S(dJ,d2),  8(dJ+1,d2+1D,  Y(dj,d2),  Y(dj  +  1,d2+1) 
be  the  corresponding  coefficients  as  defined  by  (7.13)-(7.14) . 

Using  (7. 2)- (7. 7)  and  (A.3)-(A.5)  of  Appendix  A,  we  have  for  i  »  1,2: 

Vlik+l(di+1)  "  CI  * 

X  xCVl|k(di+l5  +  Pk+1|kCdri)HT(q’k+1)R"1(q’k+-l)2Ck+l)) 

CC.l) 


Pkil|k>lCdi+1)  “  Pkil|kCdi+l)  *  Aq.k+lJR^CqJc^Hfq,^!) 


CC.2) 


And  using  (7,8)-(7.10) : 

prf-V1)  11^1. 

2 1'VTIura7 ' *"  K1)!' 


C2(k*l)-H(q,k*l)Jk+l[k(dR+1))TBCdR+1}‘1C20c>l)-H(q,k+l)5ik+1|k 


jk+1% » Tn ^ jk+1* * 1 i 


(d*+1))  -  (2Ck+l)-H(q,k+l)itk+1jk(d^i))‘BCd^‘)'i(z(k+l)- 

H(q,k*l)«krt|X*1»  (C’3) 

where  for  i  »  1,2,  8(d^+b  ■  H(q,k+l)Pk+^ j k(dk  ^)H^(q,k+l)  ♦  R(q,k+1). 
Using  (C.1)-(C.2)-(C.3) ,  after  some  tedious  calculations,  it  can 


be  shown  that: 
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APPENDIX  D 


Determination  of  the  Optimal  Gatins 


k  k 

Let  hj  and  h^  be  the  two  hypotheses  as  defined  in  Chapter  11  and 
k  k 

3(hj,h0)  be  the  corresponding  coefficient.  The  optimal  gating  is  de- 

k  k 

termined  by  zCk)  6  2  such  that  S (h^ ,hQ]  >_  0.  According  to  Chapter  11, 
we  have: 

'  *k|k-l<^  *  VhJkH£»)  -  (D.D 


^klkO'o’  *  *0 


pifk^  ■  piik.i4  ♦  hV1h 


CD.  2] 


CD. 3) 


\\W  *  P0 


CD.  4) 


PXhJ,’?)  -  PChJ’1,Yk"1)  X  Pd  X  Cl  -  P10)  x  N[z(k)  -  H^|k-1Chk), 


B(hk)3 


(D.S) 


PChJ.Y^  -  PCbJ"1,Ylt’1)  x  Pf  x  (I  -  pd)  x  Roi  x  N[zCk),I] 


CD. 6) 


sc^.hj)  -  (»k|k»^  -  ^k^Xu^  -  pk|k<’,j»'1 

.  k  .  k  PCh^Y34)  |Pk|kChi)| 

+  2  tlh+r +  *n ,  *l--t  , 
1c|k  3  *,k  0  PChJ.O  pk|kCb-)l 


CD.  7) 


Using  CD.1)-(D.6),  we  have: 
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69yhS>  ■  -  “\|k.i^)Tf-^>cxck)  -  HSkH^o-J)) 

*  -  VX  -  pklX»‘VX'k>  - 

-  HX^j ♦  zT(k)l”12Ck)  ♦  Ck)  ,  where: 

•TCk)  -  B_1Ch^  -  kJ0iJ)[Pq  -  Pk|kO»j)]"1KkChJ)  ,  CD. 8) 


PCk^”1*Yk”1) 

a .  (k)  »  2  Zn  - fa'  i  j 

J  PQi^  ,  i  j 


Pd  U-Pl0) 

+  2  Zn  -1_p  ,p  +  2  Zn  p 

U  rdjrf  *oi 


.,  l«l *  lpo' 

l»0>j)  »  pk|k0»pl 

*  '*k|k-l»j>  -  vX  -  pk|kO>i»'1(Sk|k-l(hj)  -  xo> 

CD. 9) 


Using  the  matrix  identity  CM. 3),  it  can  be  shown  that: 

Tj  00  -  CR  ♦  -  P^'V)’1  CD. 10) 

So,  BOiJ.hJ)  -  zTCk)i:'1zCk)  ♦  a .  (k) 

*  CzCk)  -  H^n.jOkJ)  -  -  PfcikCh^'1 

^Ik-lX  -  *0)T><  'X  x  t»CW  -  -  Cl*)-1. 

K^Olj) (Pq  -  ^U^Xnc-lO*^  -  x0>  *  C*k|k-1^  -  VT 
CP0  -  -  p^a^'Xii,-!^ 


Using  again  CM. 3)  it  can  be  shown  that: 
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t'0  -  t’k|k»5))'\»5’^'lKJ(,i?KPo  •  v0*)”'  * 

<po  •  pk|kchj))'1  •  lPo  -  pk|k-it,,j))‘1 


(D* 11) 


-  pkik^rl  ■  “’Mk-i^^o  •  pk|k-i°*j” 


CD. 12) 


I BOlj)  1  *  |Pk|kO>-Jl  ■  1*1  *  l\|k-i»5’I 


CD. 13) 


k.,-1 


Let  Cj(k)  -  ^0$  ♦  P^k-l^O  -  ^-1°^ 

C*k|k-l°V  '  V 

From  (D. 10) -(D. 14),  we  have: 

SOyh£)  »  CtCk)  -  HcjCk))TC-T^)C2Ck)  -  Hc.Ck)) 


CD- 14) 


♦  zTCk)r“1zCk) 


♦  2  An 


Cl-Pd)P£ 


*  2  In 


C1-P,J 


«•  An  ♦  6. Ck),  where: 


R  WJ 


POk-1,**”1) 


Si  Ck)  «  2  An - 1,  t.v- 

J  pck  \r  S 


ISjk-l^j^ 


♦  C*k|k-1°$  -  yTcPo  -  Sjk-i0^'  ^Ik-i^p  -  V 

CD. IS) 

<■>  BChJ.hS)  ■  2  2"  trir^  * 2 “  nrp  *  2"  Hf  *  *i« 

*  CzCk)  -  CT^  -  E'l)’1T^Hc.Ck))T0:‘1  -  TjH*CM  - 

CT*  -  r‘ylT*HcjC.k))  ♦  CHcJOt))T[E  -  C^)'V XCHc .  Ck))  . 


CD.  16 


Let,  E.(k)  -  TT  -  £ 

?A  a-P10)  r 

Vk)  a  2  *n  tttjtj-  +  2  *n  -t“  +  ln  fir  Vk) 

♦  (Hc.Ck))TCS  -  (T^'V^c-Ck)).  CD.  17 

Then,  SCh*  hjj)  -  c  (k)  -  UCk)  -  m. 00)TE. (k) Cz(k)  -  ®.(k)),  where 
E. (k) ,  m.Ck),  e. (k)  are  defined  by;  (D.10),  (D.14)-CD.17) . 


Q.E.D. 
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